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

更改所有复杂的类型y或类型。f尤其

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 回答