-3

I wrote this code in Fortran:

Program RK4
implicit none
real x,y,k1,k2,k3,k4,h
integer i,n
real f
read*,n,h 

x=0
y=0.001

Do i=1,n
   k1=h*f(x,y)
   k2=h*f(x+h/2.0, y+k1/2.0)
   k3=h*f(x+h/2.0, y+k2/2.0)
   k4=h*f(x+h, y+k3)
   y=y+(k1+2*k2+2*k3+k4)*(1/6.0)
   x=x+h
   write(*,*)x,y
end do
end program

real function f(x,y)
  f=1/y
end function

This code is for solving this differential equation using RK4 method:

dy/dx=1/y

I tested the code for different initial values like (x=1, y=1) and (x=0, y=0.1) and the answer was correct but for some initial values like the one in the code here (x=0,y=0.001) I do not get the right answer. Does anybody know what is the problem with the code?

Ian Bush
  • 6,996
  • 1
  • 21
  • 27
Elham Q
  • 5
  • 1
  • 3
    Why to you imagine there is a problem with the code? Depending on the problem, you will need to carefully select a step size to get a stable solution. This is likely a mathematical question – talonmies May 19 '20 at 06:33
  • 1
    At the very least you will need to tell us the incorrect values you are getting, what you think the correct values are, and what you input for n and h. – Ian Bush May 19 '20 at 07:26

1 Answers1

4

The test equation you solve, y'=1/y, has a solution y(x)=sqrt(2*x+C). For the critical initial value y(0)=0.001 you get C=1e-6, the solution has a vertical tangent at x=-5e-7, thus very close to the initial point. Such singularities will cause very large errors in numerical methods.

Or in other terms, the y-derivative of f is -1/y^2 which gives an y-Lipschitz constant of L=1e+6 for the IVP under consideration. RK4 only works at all for Lh<2.5 and sufficiently nice for Lh in the range of 1e-3 to 0.1. So the initial step size needs to be smaller than 2.5e-6, and much smaller for good results.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51