- 締切済み
常微分方程式の解き方について
常微分方程式の解き方について質問です。 Fickの第2法則の式を4次のルンゲ・クッタ法で解きたいと考えています。 4次のルンゲ・クッタ法は、 k1=f(xi,yi) k2=f(xi+h/2,yi+h/2*k1) k3=f(xi+h/2,yi+h/2*k2) k4=f(xi+h,yi+h*k3) として、 yi+1=yi+h/6*(k1+2k2+2*k3+k4)で解いています。 そこで、Fickの法則（dC/dt=D*d^2C/dx^2）を、以下の式ように変換しました。 dC/dt=D*（Cj+1+Cj）/dx^2 ←これ(dC/dt)をf(xi,yi)としています。 f(xi,yi)のyiはCとして、hはdtとして解いてみました（xiはないものとして解いています） エクセルでは、 t0 ， t0+⊿t ， ・・・ x1: Cj,n ， Cj,n+1 ， ・・・ x2: Cj+1,n ，Cj+1,n+1， ・ ・ のようなセルを組み、t0+⊿tにCn+1=Cn+h/6*(k1+2k2+2*k3+k4)の式を入れています。 上記の方法で解いてみましたが、うまくいきませんでした。 説明が分かりにくく申し訳ないですが、上記のやり方で間違っているところ、アドバスなど ありましたら回答お願いします。
- cstag
- お礼率77% (14/18)
- 物理学
- 回答数2
- ありがとう数1
- みんなの回答 （2）
- 専門家の回答
関連するQ&A
- 常微分方程式の数値解法： 陰解法のルンゲクッタ法の公式について
陽的なRunge-Kutta法は, y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4) ただし，k1=hf(xi,yi), k2=hf(xi+h/2, yi+k1/2) k3=hf(xi+h/2, yi+k2/2) k4=hf(xi+h,yi+k3) と表すことは x=[x0,xn]の範囲でルンゲ・クッタ法により数値的にy(x)を解きたいのですが， 解こうとしている問題の初期条件がx=xnの時，y=0となっており，陰的(xnからx0に向かって)に解かなければならないのだろう，と考えています． 上記の公式でhの代わりに-hを入れて，プログラムを走らせても求めたい結果に大きな差異が生じてしまい困っています． そこで，陰的なRunge-kutta法の公式には，陽的な解き方と比較してどのような修正をすればよいか，教えてください． ちなみに解きたい微分方程式は， d2y/dx2 = x と仮定します． よろしくお願いいたします！！
- 締切済み
- 数学・算数
- 連立微分方程式の解き方
d^2x/dt^2=f(t) d^2y/dt^2=g(t, y, dy/dt, z, dz/dt, d^2z/dt^2) d^2z/dt^2=h(t,d^2y/dt^2, z, dz/dt) というような連立微分方程式を解きたいのですが，どのような方法で解くことができるのでしょうか？ ルンゲ・クッタは適用できますか？ （d^2y/dt^2の式にd^2z/dt^2が変数として出てきたりしているので，わからなくなってしまいました．） 宜しくお願いします．
- ベストアンサー
- 数学・算数
- 偏微分方程式の数値解析
お忙しいところすみません． 現在，FDTD法による数値解析を行っています． FDTD法は空間及び時間微分に対して中心差分を行っていますがこれを4次のルンゲクッタ法で行うことは出来るのでしょうか？ある参考書に以下のような「連立1階微分方程式」を4次のルンゲ-クッタ法で解くというのがありました． dy/dt=f(t,y,z) dz/dt=g(t,y,z) これを4次のルンゲ-クッタ法で解くと y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4) z(n+1)=z(n)+1/6*(j1+2*j2+2*j3+j4) ここで k1=Δt*f(t(n),y(n),z(n)) j1=Δt*g(t(n),y(n),z(n)) k2=Δt*f(t(n)+Δt/2,y(n))+k1/2,z(n)+j1/2) j2=Δt*g(t(n)+Δt/2,y(n))+k1/2,z(n)+j1/2) k3=Δt*f(t(n)+Δt/2,y(n))+k2/2,z(n)+j2/2) j3=Δt*g(t(n)+Δt/2,y(n))+k2/2,z(n)+j2/2) k4=Δt*f(t(n)+Δt,y(n)+k3,z(n)+j3) j4=Δt*g(t(n)+Δt,y(n)+k3,z(n)+j3) これをマクスウェルの方程式に対応させ ε∂E(t,x,y,z)/∂t=f(t,H)=∇×H μ∂H/∂t=g(t,E)=-∇×E として出来ないこともない気がするのですが． なぜこの方法を試みようとした理由は現在，通常のFDTD定式化により数値解析を行っていますが ・プログラムは正しい ・始めの部分は正常に計算出来ている ・途中から凹凸が出始め，それが成長して発散している という事態になり，数値的不安定の可能性が強いと思ったからです．オイラー法やルンゲ-クッタ法ではこのような異常振動は進み幅を小さくすれば止まる(ということが数学的に証明されている)ようです． 何卒回答の方，よろしくお願い致します．
- 締切済み
- 物理学
- 微分方程式（ルンゲ＝クッタ）がうまくいきません
C言語で物体の自由落下運動（空気抵抗を考慮する）を微分方程式を用いて解かなければならないのですが、どうしてもグラフの形が曲線になりません。どこが間違っているのか指摘お願いします。 空気抵抗を考慮しているので、落下速度はどんどんある一定の値に近づいていくはずなのですが、グラフは直線になってしまいます。いろいろと方法を試しましたがどれも外れました。 写真は赤が今の状態（間違っています）、黒が理想形です。 ここではルンゲ＝クッタ法を使用しています。 #include <stdio.h> #include <math.h> double yd(double t,double y){ double m = 50.0; 質量 double g = 9.8; 重力加速度 double k = 0.01; 空気抵抗の定数係数 double a= g - ((k *y)/m); 運動方程式 return(a); } double runge(double t,double y,double h){ double k1,k2,k3,k4,r; double h2 = h / 2.0; k1 = yd(t,y); k2 = yd(t+h2,y+(h2*k1)); k3 = yd(t+h2,y+(h2*k2)); k4 = yd(t+h,y+(h*k3)); a = y + (h/6.0) * ((k1 + (2.0 * k2) + (2.0 * k3) + k4)); return(a); } void main(){ double y = 0.0; double h = 0.001; double t = 0.0; int i; for(i = 0; i < 1000; i++){ t += h; y = runge(t,y,h); printf("%f %f\n",t,y); } }
- ベストアンサー
- C・C++・C#
- ∑の微分について
式の導出中に∑の微分について混乱してしまいました． 考えれば考えるほど頭の中がもつれて収集がつきません． 以下のような式について考えています． 次の式の∑はどちらもiについての∑です． [ ∑{ f(xi)*log(yi) } ]dyi/dt と [ ∑{ f(xi) } * log(yj) ]dyj/dt それぞれどんな式なるのかがわかりません． 式について本当は f(xi)/yi または f(xj)/yj というものを導出したいです． そこまでが複雑すぎていろいろ理解できていないので 結果論的に考えて理解したいのですが… 混乱する中の式なので まったく意味不明な式かもしれません． 何かこれについて教えていただけないでしょうか．
- ベストアンサー
- 数学・算数
- 微分方程式のプログラムについて
次の問題を解こうとしましたが、プログラムが動かなかったので実行できませんでした。今、非常に困っています。分かりにくいプログラムですみません。 問 温度Taの空気中に置かれている物体の温度Tは、次の常微分方程式に従う。 dT/dt=-k(T-Ta) t=0でT=100[℃]、Ta=20[℃]、k=1×10^-2[s^-1]として、物体の温度が60℃にまで冷却される時間を求めよ。厳密解とも比較すること。 自分が作ったプログラムは #include <stdio.h> #include <math.h> double k = 1e-2, tmpa=20.0; double f(double t, double tmp) { return -k*(tmp-tmpa); } int main(void) { double t0 = 0.0, tmp0 = 100, t1, tmp1, ttmp1; double k1, k2, k3, k4, h; printf("刻み幅="); scanf("%lf", h); t1 = t0+h; ttmp1=tmpa+(tmp0-tmpa)+exp(-k*(t0-t0)); printf("%8.4lf %8.4lf %8.4lf\n", t0, tmp0, ttmp1); while(t0 <= 100.0){ k1 = f(t0, tmp0); k2 = f(t0+0.5*h, tmp0+0.5*h*k1); k3 = f(t0+0.5*h, tmp0+0.5*h*k2); k4 = f(t0+h, tmp0+h*k3); tmp1 = tmp0+h*(k1+2.0*k2+2.0*k3+k4)/6.0; printf("%8.4lf %8.4lf %8.4lf\n", t1, tmp1, ttmp1); if(tmp1 <= 60.0) break; t0 = t1; tmp0 = tmp1; t1 = t0+h; } return 0; } あと厳密解はT=e^-kt(T0-Ta)+Ta求められるそうです。 変な間違いが多いと思いますが、このプログラムのどこが間違っているか教えてください。厳密解だけでもいいです。お願いします。
- ベストアンサー
- C・C++・C#
- この微分方程式をお願いいたします。
スピンコートで基板上の膜厚さに関して、 添付のような計算の式が成立するときその式 （以下ではこの式を式Bと呼びます。）を用いて 微分式を使いたいですが、 私の自分の考えは、 dh/dt=式B*Δtー初期厚さ ...........式１ そして、厚さhの変化は h=初期厚さ+式１*ΔＴ と考えますが、 エクセルで換算したらなんかデータがおかしいですが、 ご存じ方がいればぜひお教授いただきませんか？ よろしくお願いいたします。 下はどんな方からもらったご回答ですが、 いま少し急いで解かせなければいけないので ご回答を頂き様のご回答を追加でもらえない状況なので、、、 また質問申し上げます。 本当に失礼しますが、 誰かがこの質問についてエクセルを用いて 関数、数式を入れて見せていただきませんか？ 本当にすみませんが、個人的に急ぐからこんなに また質問申し上げます、ぜひよろしくお願いいたします。 "添え字は略し、 f(t)=h/sqrt(1 + K*t), K=4ω^2*h^2/(3ν). とすると、 (d/dt)f(t)=(-hK/2)/sqrt( 1 + Kt)^3. ですから、導関数が必要であればこれを使ってください。 (t 以外は定数とします)" よろしくお願いいたします。
- ベストアンサー
- 数学・算数
