常微分方程式の数値解法at SIM
常微分方程式の数値解法 - 暇つぶし2ch69:名無しさん@5周年
06/07/30 15:32:42
ちょっと微分っぽく戻して

dx/dt
=f1(x, y, z)

dy/dt
=f2(x, y, z)

z=f3(x, y)
これらを差分にする

xの数列をx[n-1], x[n], x[n+1]
yの数列をy[n-1], y[n], y[n+1]
zの数列をz[n-1], z[n], z[n+1]
などとおく。
時間差を⊿tとする

前進差分なら
z[n]=f3(x[n], y[n])
x[n+1]=x[n]+⊿t*f1(x[n], y[n], z[n])
y[n+1]=y[n]+⊿t*f2(x[n], y[n], z[n])


これを順次代入し続ければ良い。これは>>65方式。

70:名無しさん@5周年
06/07/30 15:42:25
で、後退差分なら
z[n]=f3(x[n], y[n])
x[n]=x[n-1]+⊿t*f1(x[n], y[n], z[n])
y[n]=y[n-1]+⊿t*f2(x[n], y[n], z[n])


x[n]-⊿t*f1(x[n], y[n], z[n]) = x[n-1]
y[n]-⊿t*f2(x[n], y[n], z[n]) = y[n-1]

ってやるのが本来で、これを連立方程式の形にして
行列にぶち込んで解きたい所だが、

これじゃ解けないから、f1(x[n], y[n], z[n])を
テーラー展開してf1(x[n-1], y[n-1], z[n-1])と1次
差分項くらいで弄り回す必要が出てくるなあ

やっぱ後退ダメだなw
前進でも2次のホイン法あたりで解こう
最後は4次のルンゲクッタ、4次におまけがついた
ルンゲクッタジルでも使うか

71:名無しさん@5周年
06/07/31 14:50:12
オイラー法
y[n+1]=y[n]+h*f(x[n], y[n])
hが時間差で

ホイン法は、単純な常微分方程式なら

k1=h*f(x[n], y[n])
k2=h*f(x[n] + h, y[n] + k1)
y[n+1]=y[n] + (k1 + k2)/2

これを連立常微分方程式に拡張するのか

f1から出て来た値をぶち込むのがl1 , l2
f2から出て来た値をぶち込むのがm1 , m2

f3は微分方程式じゃないから動かさなくてよろしいw

z[n]=f3(x[n], y[n])

l1=h*f1(x[n], y[n], z[n])
l2=h*f1(x[n] + l1, y[n] + m1 , z[n])
x[n+1]=f1(x[n], y[n], z[n]) + (l1 + l2)/2

m1=h*f2(x[n], y[n], z[n])
m2=h*f2(x[n] + l1, y[n] + m1 , z[n])
y[n+1]=f2(x[n], y[n], z[n]) + (m1 + m2)/2

で、この問題だと時間差は1みたいだし、2回計算した平均値で
延々回してるだけじゃんかw

72:さとし
06/08/22 14:43:50
ブラジウスの公式 URLリンク(irws.eng.niigata-u.ac.jp)
F'''+(1/2)FF''=0 (境界条件:η=0 F=F'=0 および η=∞ F''=1)
をルンゲクッタで解くとあるんですけど,実際に解くときにF'',F',Fはどのように計算すれば良いでしょうか?教えて下さいm(_ _)m

73:名無しさん@5周年
06/08/22 18:51:32
>>72
一般的に
x''' + x'' + x' + x = a てな問題は
y0=x
y1=x'
y2=x''

とおくと

y0'=y1
y1'=y2
y2' + y2 + y1 + y0 = a

となって、連立の1階微分方程式となる。
1階微分方程式ならルンゲクッタで熔けるでしょ?



74:48
06/08/23 02:32:04
>>72
73氏に補足でF'(0)が未知の場合、
初期値がルンゲクッタで解くには不完全なので
射的法(シューティング法)で解く方法があります。

75:48
06/08/23 02:34:00
間違えた。F''(0)が未知の時。

76:さとし
06/08/23 21:44:25
>73、74さん
無事解決できました.ありがとうございます.m(_ _)m

77:名無しさん@5周年
07/06/16 02:20:33
てす

78:山本圭太
07/09/19 02:33:24
なんでこのスレまだ残ってるんだどーwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww

アンモナイト並の綾波だどーwwwwww

79:名無しさん@5周年
07/09/22 02:40:18
ガリレオの相対性原理は間違っていました。元日本物理学会
会長の佐藤文隆教授も昔からそれに気づき、自著に書いてお
られます。以下の文はGZKカットオフに関する記述の一部
です。

ガリレオの相対性原理は「互いに等速運動している系はみな
対等だ」というのです。
・・・・・・・・・・途中略
物理は全部、相対性原理にのっとっているといえます。
・・・・・・・・・・途中略
相対性原理でない考えは、なにか絶対系があるというもので
す。なにか特殊なものがあって、それから見てあれはあっち
に動いている、これはこっちにこう動いていると、相対運動
についてランクづけできるわけです。


ガリレオの相対性原理は、静止と等速直線運動は区別出来な
いという、慣性の法則の一部の結論によって構築されました
が、思考実験によって、400年ぶりに静止と等速直線運動は
力学的に区別出来ることが判明しました。
詳しくはこのサイトをご覧ください。
URLリンク(home9.highway.ne.jp)

80:pd3909b.osakac00.ap.so-net.ne.jp
07/10/04 14:57:58
;


81:名無しさん@5周年
08/04/09 18:09:51
「常微分方程式 予測子」でぐぐると2件目にココw

82:名無しさん@5周年
09/01/29 11:27:11
>>52
ガッ

83:名無しさん@5周年
09/05/27 22:16:16
>>36
ガッ

84:名無しさん@5周年
09/07/07 01:30:50
いや~ん!さんっていまどうしてるのだろうか

85:名無しさん@5周年
09/12/27 05:31:54
初期密度勾配 ρ(x,t=0) = C-|x| に対して
フィックの法則
dρ(x,t)/dt = Dd^2ρ(x,t)/dx^2
を使って密度ρの時間発展を観察したかった。

とりあえず差分近似でシミュってみたけど、
x=0 付近の二次導関数の急峻度が分解能によってザクザク変わって解きにくかった。
ρ(x,t=0)=C-√(1+x^2)とかちょっと解析的な関数で近似して計算するかなぁと思ったけど、
初心に戻って考えるとフィックの法則は単に
「どんどんなめらかになります」っていうだけの式だし、
x=0 付近の些細な形に囚われずに全体のおおよその形は決まって行くはず。
この初期条件ってそれなりの解析解があるのかなぁ、と思いました。
なんかいい方法教えてけれ

86:85
09/12/27 06:42:11
あ、もしかしてフーリエ変換でいけるかも。
Fourier{d/dx} = -iω より
Fourier{dρ(x,t)/dt} = Fourier{Dd^2ρ(x,t)/dx^2} = -Dω^2 Fourier{ρ(x,t)}
Fourier{ρ(x,t)} = exp(-Dtω^2)Fourier{ρ(x,0)}
フーリエ逆変換して
ρ(x,t) = (2Dt)^(1/2)exp(-x^2/(4Dt))*ρ(x,0)
(*は x に関する畳み込み演算)

うわああああ。時間比例する分散をもつガウス分布の畳み込みになるんだ。
これ初期配置どんな形でも解けるじゃん。
っていうかなんか中心極限定理そのままじゃん。
おれはいまはげしく感動しているぞ…

というか、ここまで書けよ、フィックの野郎…

87:85
09/12/27 06:52:52
逆フーリエ変換ちょっと間違ってたな。
ρ(x,t) = (4πDt)^(1/2)exp(-x^2/(4Dt))*ρ(x,0)

つーかフィックの法則・正規分布でぐぐったらいっぱい出てくるわ…
まあ自分で解いて感動できたのでOK

88:85
09/12/27 15:56:48
さらに訂正
ρ(x,t) = (4πDt)^(-1/2)exp(-x^2/(4Dt))*ρ(x,0)

89:名無しさん@5周年
10/04/17 12:20:43
hoshu


最新レス表示
レスジャンプ
類似スレ一覧
スレッドの検索
話題のニュース
おまかせリスト
オプション
しおりを挟む
スレッドに書込
スレッドの一覧
暇つぶし2ch