xsr
xsr

Reputation: 13

Fortran array unclassiable statement

  program rk4
  real x(200), y(200), h, k1, k2, k3, k4
  integer i
  read*,x(0),y(0),h
  do i=0,10
  k1=x(i)-y(i)+2
  x(i+1)=x(i)+h
  k2=(x(i)+h/2)-(y(i)+(h/2)*k1))+2
  k3=(x(i)+h/2)-(y(i)+(h/2)*k2))+2
  k4=(x(i+1))-(y(i)+h*(k3))+2
  y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4)
  print*, 'x=', x(i+1), 'y=', y(i+1)
  enddo
  end

At line 9 and 10:

  k2=(x(i)+h/2)-(y(i)+(h/2)*k1))+2
  k3=(x(i)+h/2)-(y(i)+(h/2)*k2))+2

I'm getting "Unclassifiable statement at (1)", with (1) pointing to k2 and k3. I can't see what I'm doing wrong as k1 and k4 follows a similar structure and it seems like there isn't anything wrong with them.

Upvotes: 1

Views: 87

Answers (1)

roygvib
roygvib

Reputation: 7395

It seems that the error message comes from too many right parentheses ")" for k2 and k3. Another error is that arrays x and y need to be declared as x(0:200) and y(0:200), because you are accessing x(0) and y(0). If the above two points are fixed, the code should work correctly.

As side notes, I really recommend to put implicit none, which is useful to detect potential bugs, and use floating-point literals such as 2.0 rather than 2 in floating point arithmetic (except for powers like x**2). In the code below, I have compared the analytical solution with your RK4 result, which seem to agree well with each other.

program rk4
implicit none                                   !<--- this is useful
real x(0:200), y(0:200), h, k1, k2, k3, k4      !<--- indices now start from 0
integer i
read *, x(0), y(0), h

do i = 0, 10
    x(i+1) = x(i) + h
    k1 = x(i) - y(i) + 2.0                                  ! y' = x - y + 2
    k2 = ( x(i) + h/2.0 ) - ( y(i) + h/2.0 * k1 ) + 2.0     !<--- ")" has been removed
    k3 = ( x(i) + h/2.0 ) - ( y(i) + h/2.0 * k2 ) + 2.0     !<--- here also
    k4 = ( x(i+1)       ) - ( y(i) + h     * k3 ) + 2.0
    y(i+1) = y(i) + h/6.0 * ( k1 + 2.0*k2 + 2.0*k3 + k4 )
    print*, 'x=', x(i+1), 'y(rk4)=', y(i+1), &
            'y(ana)=', x(i+1)+1.0 + (y(0)-x(0)-1.0) * exp(-x(i+1)+x(0))
enddo
end

Upvotes: 1

Related Questions