Solving first order complex ODE

Question

I have a question about solving the first order complex differential equation. I used Runge-Kutta but the answer does not seem to be right.

This is my equation:

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

The results for ODEs are OK but for complex equations it does not seem to be correct.

I compiled my code with gfortran and I plotted the diagram with gnuplot for ordinary differential equation without i. The diagram was correct for the solution of the equation. but I added i in my equation and the diagram was just a straight line. Also, gnuplot just draw the real part.

The code was compiled without any compile-time errors.

This is my code:

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

and the output:

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

If you have a problem with Runge-Kutta, post a question describing the details of your problem and your code. — Vladimir F Героям слава, Jan 20 '16 at 13:58
That link to your code is certainly not enough. The code must be inside the question. You must describe the problem in detail. What are your results? Which results do you expect? Are there any error messages? How do you compile the code? — Vladimir F Героям слава, Jan 20 '16 at 14:13
Imaginary unit is (0,1). You should readily see the error in `f4=h*f(t+h,x+f3)`. — Lutz Lehmann, Jan 20 '16 at 14:28
Obviously, if you want to compute with complex numbers, this should be represented in the data type. All but the time and time step variables should have type `complex`. You should have gotten a compiler warning about incompatible types at assignment. — Lutz Lehmann, Jan 20 '16 at 15:08
Do you mean I have to change all of the variables such as a, b etc. except of time and time step? — shimool, Jan 20 '16 at 15:18
a and b are the endpoints of the time interval, thus fall under time. And yes, everything else that was `real` should be `complex`. — Lutz Lehmann, Jan 20 '16 at 15:21
In your function `f` you have a local integer variable `i` which is never assigned to or initialized but appears somewhat important to the result. Although you use `implicit none` in the main program you don't in the other scopes. — francescalus, Jan 20 '16 at 16:02
If learning to program, and hopefully also later, always compile with option `-Wall` to display all warnings. In such a short program there should be no warnings. -- It is doubtful that the real version worked correctly since the y-step is always from the initial fixed value `y1` and not the current `y`. — Lutz Lehmann, Jan 20 '16 at 16:02

Lutz Lehmann · Answer 1 · 2016-01-20T16:36:53.707

Change everything that is y or f to type complex. Especially

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

and correct the RK4 loop to actually update the current value of y (and open the output file only once before the loop)

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

By the way, the exact solution can be computed as

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

Solving first order complex ODE

1 Answers1