1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Fortran programming to solve linear equation for ode

  1. Dec 1, 2013 #1
    Fint the exact solution of the system
    dy/dt = -15y-25z
    dz/dt=-47y-85z
    with inital condition y(0)=2, z(0)=5
    either by writing the equation in matrix form as dx/dt =AX where x=(y z) and diagonalising the matrix A, or otherwise.
    Using fortran programming with second order adam bashforth method, plot the numerical soltuions for y_i and z_i against the exact solution between t=0 and t=0.1 for different time steps h(large and small), chosen to distinguish between instability and stability(or in the case of an unconditionally stable scheme, between being oscillatory and non oscillatory)
    Here is my code

    !Second order Adams-Bashforth method for ODE
    !dy/dt= -15y-25z
    !dz/dt=-47y-85z
    ! with intial condition y(0)=2, z(0)=5

    Program adamstwo
    Implicit None
    Real, allocatable :: y(:),t(:), z(:), texact(:),yexact(:),zexact(:)
    Real:: tend, h, k1,k2,k3,k4
    Real, parameter :: exp=2.718282
    Integer:: NI, i
    Print*, 'Enter the final time '
    read*, Tend
    Print*, 'Enter number of timesteps to take'
    read*, NI
    h= Tend/NI
    Print*, 'This gives stepsize h=',h
    allocate (t(0:NI), y(0:NI),z(0:NI))
    !Initial Conditions
    t(0) = 0
    y(0) = 2
    z(0) = 5

    !After using runge kutta method, I found out k1 =-155 and k2= 3860*h-155 ,
    k1 = -155
    k2 = 3860*h-155
    !we know that y(n+1) =y(n) + h/2(k1+k2)at n=0
    y(1) = y(0) + h/2 *( k1 + k2)

    !After using runge kutta method, I found out k3 =-155 and k4= 3860*h-155
    k3= -519
    k4= 44068*h-519
    !we know that z(n+1) =z(n) + h/2(k3+k4) at n=0
    z(1) = z(0) + h/2 *( k3 + k4)

    t(1) = h

    open(10,file='adams_two results.m')
    ! Loop through the number of steps to calculate the following at each step
    do i = 2, NI
    t(i) = i*h

    !Second order Adam bashforth for all n
    y(i+1) = y(i)+ (h/2)*(3*(-15*y(i)-25*z(i))-(-15*y*(i-1)-25*z*(i-1)))
    z(i+1)= z(i)+ (h/2)*(3*(-47*y(i)-85*z(i))-(-47*y*(i-1)-85*z*(i-1)))

    end do

    !Print out the Approximate solution
    write(10,*) 'ApproximateSolution =[', t(0),y(0),z(0)
    do i = 0, NI
    write(10,*) t(i),y(i),z(i)
    end do
    write(10,*) t(NI), y(NI),z(NI),']'
    allocate (texact(0:NI), yexact(0:NI),zexact(0:NI))
    texact(0)=0
    yexact(0)=2
    zexact(0)=2
    do i = 1, NI
    texact(i) = i*h
    yexact(i)=1/(4*(i*h)**4 +1)
    yexact(i)= 0.4385* exp**((-50+20*sqrt(6.0))*(i*h)) + 1.5613 *exp**((-50-20*sqrt(6.0))*(i*h))
    zexact(i)= -0.2455* exp**((-50+20*sqrt(6.0))*(i*h)) + 5.2453* exp**((-50-20*sqrt(6.0))*(i*h))
    end do

    !Print out the exact solution
    write(10,*) 'ExactSolution = [',texact(0), yexact(0),zexact(0)
    do i = 0, NI
    write(10,*) texact(i), yexact(i),zexact(i)
    end do
    write(10,*) texact(NI), yexact(NI),zexact(NI),']'
    write(10,*) "plot(ApproximateSolution(:,1),ApproximateSolution(:,2),'g',ExactSolution(:,1),ExactSolution(:,2),'r')"
    write(10,*) "xlabel('time'),ylabel('y'),legend('Approximate AB[2] Solution','Exact Solution')"
    close(10)

    end program adamstwo

    I got the run time error :
    Error 112,Reference to undefined variable,array element or function result(/UNDEF)
    main - in file adamstwo.f95 at line 44 [+0877] ---
    (i don't have matlab file because of this runtime in fortran programming.)

    Please help.
    (Note:I have calculated the exact solution of the system by hand which is correct)
     
  2. jcsd
  3. Dec 1, 2013 #2

    Mark44

    Staff: Mentor

    Here are lines 44 and 45 of your code. I have included both lines, because even if you fix line 44, the compiler will still give you an error for line 45
    Code (Text):
    y(i+1) = y(i)+ (h/2)*(3*(-15*y(i)-25*z(i))-(-15*[color="red"]y*(i-1)[/color]-25*[color="red"]z*(i-1)[/color]))
     z(i+1)= z(i)+ (h/2)*(3*(-47*y(i)-85*z(i))-(-47*[color="red"]y*(i-1)[/color]-85*[color="red"]z*(i-1)[/color]))
     
    There are no y and z variables in your code - there are arrays named y and z. Where you have y * <something> and z * <something> in the lines above, you should have y(i - 1) and z(i - 1), to access the element in each array at index i - 1.

    It's possible that there are other errors, but if the fix the above two lines, that should take care of the run-time error you're seeing.
     
  4. Dec 2, 2013 #3
    ok i changed it
    y(i+1) = y(i)+ (h/2)*(3*(-15*y(i)-25*z(i))-(-15*y(i-1)-25*z(i-1)))
    z(i+1)= z(i)+ (h/2)*(3*(-47*y(i)-85*z(i))-(-47*y(i-1)-85*z(i-1)))
    and added this code : t(NI)=NI *h and its working. Thank you mark 44 and others too.
     
  5. Dec 2, 2013 #4
    i want to plot two graph : yn against exact solution and zn against exact solution separately
    Here is fortran code to open file for matlab:
    write(10,*) "plot(yn(:,1),yn(:,2),'b',ExactSolution(:,1),ExactSolution(:,2),'r')"
    write(10,*) "xlabel('time'),ylabel('y'),legend('yn','Exact Solution')"
    write(10,*) "plot(zn(:,1),zn(:,2),'g',ExactSolution(:,1),ExactSolution(:,2),'r')"
    write(10,*) "xlabel('time'),ylabel('y'),legend('zn','Exact Solution')"
    But it does not work when i open the matlab file , it gives the numerical solution of yn, zn and exact solution. but when i run to plot the two graphs , it only plot one graph which is zn aginst exact solution. it does not open or plot yn against exact solution.
     
    Last edited: Dec 2, 2013
  6. Dec 2, 2013 #5
    !Second order Adams-Bashforth method for ODE
    !dy/dt= -15y-25z
    !dz/dt=-47y-85z
    ! with intial condition y(0)=2, z(0)=5
    !Rabindra Gurung Qs part2
    Program adamstwo
    Implicit None
    Real, allocatable :: y(:),t(:), z(:), texact(:),yexact(:),zexact(:)
    Real:: tend, h, k1,k2,k3,k4
    Real, parameter :: exp=2.718282
    Integer:: NI, i
    Print*, 'Enter the final time '
    read*, Tend
    Print*, 'Enter number of timesteps to take'
    read*, NI
    h= Tend/NI
    Print*, 'This gives stepsize h=',h
    allocate (t(0:NI), y(0:NI),z(0:NI))
    !Initial Conditions
    t(0) = 0.0
    y(0) = 2.0
    z(0) = 5.0

    !After using runge kutta method, I found out k1 =-155 and k2= 3860*h-155 ,
    k1 = -155.0
    k2 = 3860.0*h-155.0
    !we know that y(n+1) =y(n) + h/2(k1+k2)at n=0
    y(1) = y(0) + h/2 *( k1 + k2)

    !After using runge kutta method, I found out k3 =-155 and k4= 3860*h-155
    k3= -519.0
    k4= 44068.0*h-519.0
    !we know that z(n+1) =z(n) + h/2(k3+k4) at n=0
    z(1) = z(0) + h/2.0 *( k3 + k4)

    open(10,file='adams_two results.m')
    ! Loop through the number of steps to calculate the following at each step
    do i = 1, NI-1
    t(i) = i*h

    !Second order Adam bashforth for all n
    y(i+1) = y(i)+ (h/2.0)*(3.0*(-15.0*y(i)-25.0*z(i))-(-15.0*y(i-1)-25.0*z(i-1)))
    z(i+1)= z(i)+ (h/2.0)*(3.0*(-47.0*y(i)-85.0*z(i))-(-47.0*y(i-1)-85.0*z(i-1)))
    end do
    t(NI)=NI *h

    !Print out the yn
    write(10,*) 'yn =[', t(0),y(0)
    do i = 0, NI
    write(10,*) t(i),y(i)
    end do
    write(10,*) t(NI), y(NI),']'

    !Print out the zn
    write(10,*) 'zn =[', t(0),z(0)
    do i = 0, NI
    write(10,*) t(i),z(i)
    end do
    write(10,*) t(NI), z(NI),']'

    allocate (texact(0:NI), yexact(0:NI),zexact(0:NI))
    texact(0)=0
    yexact(0)=2
    zexact(0)=2
    do i = 1, NI
    texact(i) = i*h
    yexact(i)= 0.4385* exp**((-50.0+20.0*sqrt(6.0))*(i*h)) + 1.5613 *exp**((-50.0-20.0*sqrt(6.0))*(i*h))
    zexact(i)= -0.2455* exp**((-50.0+20.0*sqrt(6.0))*(i*h)) + 5.2453*exp**((-50.0-20.0*sqrt(6.0))*(i*h))
    end do
    !Print out the exact solution
    write(10,*) 'ExactSolution = [',texact(0), yexact(0),zexact(0)
    do i = 0, NI
    write(10,*) texact(i), yexact(i),zexact(i)
    end do
    write(10,*) texact(NI), yexact(NI),zexact(NI),']'
    write(10,*) "plot(yn(:,1),yn(:,2),'b',ExactSolution(:,1),ExactSolution(:,2),'r')"
    write(10,*) "xlabel('time'),ylabel('y'),legend('yn','Exact Solution')"
    write(10,*) "plot(zn(:,1),zn(:,2),'g',ExactSolution(:,1),ExactSolution(:,2),'r')"
    write(10,*) "xlabel('time'),ylabel('y'),legend('zn','Exact Solution')"
    write(10,*)"hold all"
    close(10)
    end program adamstwo

    (Here is my final corrected fortran code and I would like to plot yn against exact solution and zn against exact solution separately but it giving only one plot in matlab.)
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted