2

一階複素微分方程式の解法について質問があります。Runge-Kutta を使用しましたが、答えが正しくないようです。

これは私の方程式です：

y'=exp(-2*t)-i*y

ODE の結果は問題ありませんが、複雑な方程式の場合は正しくないようです。

コードを gfortran でコンパイルし、gnuplot を使用して常微分方程式を使用せずに図をプロットしましたi。方程式の解について、図は正しかった。しかし、方程式を追加iすると、図は単なる直線になりました。また、gnuplot は実部を描画するだけです。

コードは、コンパイル時エラーなしでコンパイルされました。

これは私のコードです:

program kutta
implicit none
real:: a=0.0, b=8.0,y=0.1,h,t
integer::n=40
complex::i=(0,1)


interface



function f(t,y)
real::t,y
end function f
end interface

h=(b-a)/real(n)
t=a
call rk4(f,t,y,h,n)
end program kutta



function f(t,y)
f=(exp(-2*t)-2*y*i)
end function f

subroutine rk4(f,t,y,h,n)
real::f1,f2,f3,f4,t1,y1

real::t,y,h

integer::k

interface
function f(t,y)
real::t,y
end function f
end interface




t1=t
y1=y
do k=1,n
f1=h*f(t,y)
f2=h*f(t+h/2.0,y+f1/2)
f3=h*f(t+h/2.0,y+f2/2)
f4=h*f(t+h,y+f3)
y=y1+(f1+2*f2+2*f3+f4)/6.0
t=t1+h*real(k)

open(unit=1,file='data.dat')
write(1,*) t, y
!print*,k,t,y
end do
close (unit=1)


end subroutine rk4

そして出力：

0.200000003       1.22892008E+14
  0.400000006       1.46381659E+29
  0.600000024                  NaN
  0.800000012                  NaN
   1.00000000                  NaN
  ....
4

1 に答える 1

2

yorfを複雑なタイプに変更します。特に

complex function f(t,y)
  real::t
  complex::y
  complex::i=(0,1)
  f = exp(-2*t)-i*y 
end function f

RK4ループを修正して、現在の値を実際に更新しますy（そして、ループの前に一度だけ出力ファイルを開きます）

do k=1,n
  f1=h*f(t      ,y       )
  f2=h*f(t+0.5*h,y+0.5*f1)
  f3=h*f(t+0.5*h,y+0.5*f2)
  f4=h*f(t+    h,y+    f3)

  y=y+(f1+2*f2+2*f3+f4)/6.0
  t=t+h

  write(1,*) t, y
  !print*,k,t,y
end do

ちなみに、正確な解は次のように計算できます。

complex function yexact(t,y0)
  real::t
  complex::y0
  complex::i=(0,1)
  complex::A=(0.4,0.2)
  yexact = -A*exp(-2*t) + (A+y0)*exp(-i*t)
end function yexact
于 2016-01-20T16:07:55.047 に答える