Euler Method for Numerical Solutions

From a Taylor expansion, we find the LOCAL error is . The order of a numerical method for solving ODEs is the Global error - for Euler (i.e., first order accuracy).

To get better accuracy, one could simply include more terms in the Taylor series expansion. But this approach, although providing better formal order-of-accuracy, is cumbersome. We will study more flexible approaches.

Here is a simple code for implementing the Euler algorithm.

         program odesolve

      Implicit none
      Double precision y(1001),t(1001)
      Double precision deltat
      Integer i,nsteps

      nsteps=100
      t(1)=0.d0
      deltat=0.01d0
      y(1)=1.d0


      call euler(t,y,deltat,nsteps)


      do i=1,nsteps
c        write(*,*) t(i),y(i)
        write(1,1412) t(i),y(i)
      enddo
1412  format('',f9.6,1x,f12.7)


      stop
      end

c     *******************************************************
      subroutine euler(t,y,deltat,nsteps)
      Implicit none
      Double precision y(1001),t(1001)
      Double precision deltat,ff
      Integer i,nsteps
      External ff

      do i=2,nsteps
        y(i)=y(i-1)+deltat*ff(t(i-1),y(i-1))
        t(i)=t(i-1)+deltat
      enddo

      return
      end
c     *******************************************************
      double precision function ff(s,z)
      Implicit none
      Double precision s,z

      ff=-z

      return
      end

Run the code above, using different functions, and plot the solutions using the tool xmgr (on Unix machines); an older version is called xvgr. The code should produce a file `fort.1'. The command `xmgr fort.1' should open a window for you and plot your solution. Play with the xmgr window - the plot menu lets you put titles and legends, change symbols, etc., and the file menu will let you save the file (or a postscript version of the plot), print it out, and a couple of other things. Get comfortable with it - we will be using this as a basis for studying other algorithms.