• 締切済み

偏微分方程式の数値解析

お忙しいところすみません. 現在,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定式化により数値解析を行っていますが ・プログラムは正しい ・始めの部分は正常に計算出来ている ・途中から凹凸が出始め,それが成長して発散している という事態になり,数値的不安定の可能性が強いと思ったからです.オイラー法やルンゲ-クッタ法ではこのような異常振動は進み幅を小さくすれば止まる(ということが数学的に証明されている)ようです. 何卒回答の方,よろしくお願い致します.

みんなの回答

  • shkwta
  • ベストアンサー率52% (966/1825)
回答No.1

刻み幅を小さくしても発散するのでしょうか?(光速との関係で、刻み幅が大きいと本質的に計算不可能なため)

daitasuki
質問者

お礼

回答ありがとうございます. FDTDにはCourant(クーラン)条件がありますが その条件の50分の1という風にかなり小さくしても 発散が起きてしまいます.

関連するQ&A

  • 常微分方程式の解き方について

    常微分方程式の解き方について質問です。 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)の式を入れています。 上記の方法で解いてみましたが、うまくいきませんでした。 説明が分かりにくく申し訳ないですが、上記のやり方で間違っているところ、アドバスなど ありましたら回答お願いします。

  • 微分方程式(ルンゲ=クッタ)がうまくいきません

    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); } }

  • 常微分方程式の数値解法: 陰解法のルンゲクッタ法の公式について

    陽的な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が変数として出てきたりしているので,わからなくなってしまいました.) 宜しくお願いします.

  • ルンゲクッタ法の二階微分方程式(Fortran)

    数値計算の演習問題で以下の二階微分方程式をルンゲ・クッタで解けという問題があります。 -y"+x^2・y=e・y(eは定数、”・”は単なる掛け算) y(0)=1, y'(0)=0, 0<=x<=2までを計算せよ。 これは y’=z・・・(1)   z'=(x^2-e)y・・・(2) この2つの連立方程式を解けばよいところまではわかります。 まず(2)を解くときにルンゲ・クッタの場合 (k1+2k2+2k3+k4)/6の項(←公式の右辺第二項)のk(1~4)を求めなければいけません。 質問はkの求め方です。 本にはy'=f(x,y,z) , z'=g(x,y,z)とおけば (2)の場合だと例えばk1は k1=g(xn,yn,zn)dxで計算する。と書いてあります しかしz'=(x^2-e)y(←zが入ってない) なので、計算すると k1=g(xn,yn)dxとなってしまうんですがどうなんでしょう? おそらくどこかで勘違いしてると思うんです。 長い質問になってしまいましたがどうかご教授のほどよろしくお願いします。

  • 微分方程式の逆算

    以下の微分方程式について、Excelで計算しようと思っています。 dy/dt = A f(t) - g(t) 普通はこのような問題の場合、 A、f(t)、g(t)を定め、ルンゲ=クッタ法などを用いてyの終端値を求めると思います。 しかし今回求めたいのは、yではなくAの値です。 f(t)、g(t)、yの終端値が決まっており、Aの値を逆算したいのです。 とにかくしらみ潰しに値を入れれば最適値が見つけられるのでしょうが それでは時間がかかりすぎるのではないかと懸念しています。 どのような解法があるか、教えてください。

  • 常微分方程式の問題で。

    常微分方程式の問題で オイラー法、ルンゲ・クッタ法、ルンゲ・クッタ・ギル法、ミルン法、アダムス法より2つの方法を用いて解け。 y''+y'/x+y=0 初期条件y'/x=-1/2 0<=x<=4 h=0.4 という問題なのですが、基本的に解き方がわかりません。 二階微分なので連立ルンゲ・クッタで解くのかな?っと思ったのですが 結局それからどうすればいいのかわかりません。 2つの方法とありますが、これはどの方法でも解けるものなのでしょうか? 本当はC言語のプログラムを組まなければいけないのですが、 この理論がわからないのでおねがいします。

  • 数値解析法

    このHeun法のプログラムをRunge-Kutta法にするにはどうしたらいいですか? #include <stdio.h> #include <math.h> double f1(double y) { return y; } double f2(double y) { return -4*y; } int main(){ double a=0; double b; int m=10; int n; double h; double x,y; int k; double e; double f; double k1,k2; printf("Heun法計算例:y=e^x, y=1/e^4x\n\n"); // y = e^x b = 1; for(n=100;n<=10000;n*=100){ h = (b-a)/n; printf("y' = y: h(=dx) = %.1e (y=e^x)\n",h); x = a; y = 1; for(k=0;k<=n;k++) { x = k*h; if(k%(n/m)==0) { f = exp(x); e = fabs(y-f); printf("x=%.2f, y=%f, e^x=%f er=%.0e\n",x,y,f,e); } // Heun's method k1 = h*f1(y); k2 = h*f1(y+k1); y += (k1+k2)/2; } }

  • 常微分方程式

    一階の常微分方程式の解を求めるプログラムを書いてみたのですが、分からない事柄があったので質問させていただきました。 なお、このプログラムは dy/dx = e^(3x)+2y という微分方程式の解をオイラー法とホイン法で計算する(厳密値も計算している)というものです。 このプログラムを実行してみたところ、オイラー法で計算した値は正しいことが分かりました。しかし、ホイン法で計算した値が微妙に違っているようでした。(出力結果の模範解答と照らし合わせた結果) 個人的にはプログラムの内容は正しいと思ったのですが、どこかに間違いがあるのでしょうか? どなたか分かる方がいらっしゃいましたら、教えていただけると幸いです。 Option Explicit Sub main() Dim h As Double, xmax As Double Dim x As Double, y As Double Dim yS As Double, yE As Double, yH As Double Dim k1 As Double, k2 As Double Dim f As Double, g As Double Dim n As Double x = 0 yS = 1 yE = 1 yH = 1 h = 0.01 n = 2 xmax = 1 Sheet1.Cells(2, 1) = x Sheet1.Cells(2, 2) = yE Sheet1.Cells(2, 3) = yH Sheet1.Cells(2, 4) = yS With Worksheets("Sheet1") Do yS = Exp(3 * (x + h)) rhs x, yE, f yE = yE + h * f rhs x, yH, f rhs2 x, yH, h, k1, g k1 = h * f k2 = h * g yH = yH + (k1 + k2) / 2 x = x + h n = n + 1 .Cells(n, 1) = x .Cells(n, 2) = yE .Cells(n, 3) = yH .Cells(n, 4) = yS Loop Until x >= xmax End With End Sub Sub rhs(x As Double, y As Double, f As Double) f = Exp(3 * x) + 2 * y End Sub Sub rhs2(x As Double, y As Double, h As Double, k1 As Double, g As Double) g = Exp(3 * (x + h)) + 2 * (y + k1) End Sub Sub clear() Sheet1.Columns(1).ClearContents Sheet1.Columns(2).ClearContents Sheet1.Columns(3).ClearContents Sheet1.Columns(4).ClearContents Sheet1.Cells(1, 1) = "xの値" Sheet1.Cells(1, 2) = "オイラー法" Sheet1.Cells(1, 3) = "ホイン法" Sheet1.Cells(1, 4) = "厳密解" End Sub

  • ルンゲクッタ法(RK4)で斜方投射の問題を解く

    下のように自由落下する物体の速度の変化を求めるプログラムをC言語で作成しました。これを応用して斜方投射された物体のx、y座標の「位置」の変化を求めるプログラムをルンゲクッタ法(RK4)で作成したいです。x,yについて2つの式を立てないといけないと思うのですがどのような式を立ててプログラムを組めかよいかよく分かりません。どなたか式とできればプログラムを教えて下さい。 k:空気抵抗、m:物体の質量 g:重力加速度としています。 #include <stdio.h> double yd(double t,double y){ double k=0.1; double m=0.1; double g=9.8; double r=g-((k*y)/m); return r; } 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)); r=y+(h/6.0)*((k1+(2.0*k2)+(2.0*k3)+k4)); return r; } void main(){ double y=0.0; double h=0.001; double t=0.0; int i; for(i=0;i<100;i++){ t+=h; y=runge(t,y,h); printf("%f %f\n",t,y); } }