{"id":183,"date":"2017-09-13T10:55:07","date_gmt":"2017-09-13T13:55:07","guid":{"rendered":"http:\/\/www.professores.uff.br\/diomarcesarlobao\/?page_id=183"},"modified":"2017-09-13T10:55:07","modified_gmt":"2017-09-13T13:55:07","slug":"odeint-f","status":"publish","type":"page","link":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/odeint-f\/","title":{"rendered":"odeint.f"},"content":{"rendered":"
c<html>\r\nc<head><title>odeint.f<\/title><\/head>\r\nc<body>\r\nc<pre>\r\nc<a name=\"module\"><font color=\"FF0900\">\r\n      module precision\r\nc<\/font><\/a>\r\nc<a name=2><font color=FF0000>\r\n         integer, parameter :: rp=selected_real_kind(13,300)\r\nc<\/font>\r\n      end module precision\r\n      program odeint\r\n      use precision\r\n      implicit none\r\nc\r\nc  Demonstration of Runge-Kutta, Adams-Bashford, and Adams-Moulton for\r\nc  solution of an Ordinary Differential Equation in the form:\r\nc\r\nc     d y\r\nc     ---  =  f ( x , y )\r\nc     d x\r\nc\r\nc  The arrays below allow storage of the results for the four most recent\r\nc  time steps.  For example y(0) contains the value of y at the beginning\r\nc  of the current time step, y(1) contains the value of y at the beginning\r\nc  or the previous time step, y(2) contains the value of y 2 steps back, etc.\r\nc\r\n      real (rp) f(0:3),y(0:3),fa(0:3),ya(0:3),fimp(0:3),yt(0:3),ft(0:3)\r\n      real (rp) h, ynew, yrk, yab, yanew, xnew,  dfdy, xn, xk1,\r\n     &          xk2, xk3, xk4, yam, yimpn, yimp, gy, dgdy, dely, yabt,\r\n     &          fnew, func, ftnew\r\n      integer nsteps, istep, i, it\r\n      character*12 filename\r\n      character*5 xstep\r\nc   Read in the integration step size\r\nc <a name=\"write\"><font color=\"FF0000\">\r\n    1 write(*,'(a)',advance='no') 'Provide Integration Step Size: '\r\nc <\/font><\/a> \r\n      read (*,*)h\r\nc   Make up output a file name that means something\r\nc   Combining 'odeout-' with the step size\r\n      write(xstep,'(f5.3)') h\r\n      filename = 'odeout-'\/\/adjustl(xstep)\r\n      open(10,file=filename)\r\n      write(10,2010)\r\n 2010 format(4x,'x',10x,'correct',5x,'Runge',7x,'Adams',7x,'Adams',4x,\r\n     $  'Implicit',\/15x,'answer',6x,'Kutta',7x,'Bashford',4x,'Moulton',\r\n     $  2x,'Adams-Moultn',\/)\r\nc\r\nc     Get A number of integration steps\r\nc\r\n      write(*,'(a)',advance='no') 'Number of integration steps: '\r\n      read(*,*) nsteps\r\nc\r\nc    Take at least 5 steps\r\nc\r\n      nsteps=max(nsteps,5)\r\nc\r\nc    Here We hardwire the initial conditions\r\nc\r\n      ynew=0.\r\n      xnew=0.\r\n      call answer(xnew,yanew)\r\nc\r\nc   Initialize values of previous timesteps\r\nc\r\n      do 10 i=1,3\r\n      y(i)=0.\r\n      ya(i)=0.\r\n      yt(i)=0.\r\n      f(i)=0.\r\n      fimp(i)=0.\r\n      ft(i)=0.\r\n   10 fa(i)=0.\r\nc\r\nc   Do Three steps of Runge Kutta\r\nc\r\n      do 40 istep=1,3\r\nc\r\nc   Shuffle past data\r\nc\r\n      \tdo 30 i=3,1,-1\r\n           yt(i)=yt(i-1)\r\n           ya(i)=ya(i-1)\r\n      \t   y(i)=y(i-1)\r\n           f(i)=f(i-1)\r\n           ft(i)=ft(i-1)\r\n           fimp(i)=fimp(i-1)\r\n   30      continue\r\n        y(0)=ynew\r\n        yt(0)=ynew\r\n        ya(0)=yanew\r\n      \txn= xnew\r\nc\r\nc    Make evaluations for this time step\r\nc\r\n        f(0)=func(xn,y(0),dfdy)\r\n        ft(0)=f(0)\r\nc\r\nc   Initialize numbers needed for the Implicit Adams-Moulton method.\r\nc   Note that I am cheating here  fimp should be\r\nc   loaded with the results of an implicit Runge-Kutta\r\nc   or explicit formulation for variable step size using\r\nc   smaller steps for initialization\r\nc\r\n        fimp(0)=func(xn,ya(0),dfdy)\r\nc\r\nc   Actual  Runge-Kutta step\r\nc\r\n      \txk1=h*func(xn,y(0),dfdy)\r\n      \txk2=h*func(xn+.5*h,y(0)+.5*xk1,dfdy)\r\n      \txk3=h*func(xn+.5*h,y(0)+.5*xk2,dfdy)\r\n      \txk4=h*func(xn+h,y(0)+xk3,dfdy)\r\nc\r\nc   Values at the new time level (end of time step)\r\nc   If you look at this carefully it is a close relative of a\r\nc   Simpson's Rule integration\r\nc\r\n        ynew=y(0)+(xk1+2.*xk2+2.*xk3+xk4)\/6.\r\n        xnew=xn+h\r\nc\r\nc   Get the actual answer for comparison\r\nc\r\n        call answer(xnew,yanew)\r\n        write(10,2000)xnew,yanew,ynew\r\n   40 continue\r\n 2000 format(1p,6e12.5)\r\nc\r\nc    Now that we have starter values for Adams-Bashford and\r\nc    Adams-Moulton, Finish integration with all methods\r\nc\r\nc    Current y for Runge-Kutta\r\nc\r\n      yrk=ynew\r\nc\r\nc    Current y for Adams-Bashford\r\nc\r\n      yab=ynew\r\nc\r\nc    Current y for Adams-Moulton\r\nc\r\n      yam=ynew\r\nc\r\nc    Current y for Implicit Adams-Moulton\r\nc\r\n      yimpn=yanew\r\nc\r\nc\r\nc\r\n      do 80 istep=4,nsteps\r\nc\r\nc   Set up for next time step by shifting stored values of previous time\r\nc   steps to keep only the results of the three previous steps\r\nc\r\n      \tdo 60 i=3,1,-1\r\n           ya(i)=ya(i-1)\r\n           yt(i)=yt(i-1)\r\n      \t   y(i)=y(i-1)\r\n           f(i)=f(i-1)\r\n           ft(i)=ft(i-1)\r\n           fimp(i)=fimp(i-1)\r\n   60      continue\r\n      \txn= xnew\r\n        y(0)=yam\r\n        yt(0)=yab\r\n        f(0)=func(xn,y(0),dfdy)\r\n        ft(0)=func(xn,yt(0),dfdy)\r\nc\r\nc   The Implicit Adams-Moulton requires a Newton Iteration (DO loop 70)\r\nc   to solve a non-linear equation in the new-time value of y\r\nc\r\n        yimp=yimpn\r\n        fimp(0)=func(xn,yimp,dfdy)\r\n        xnew=xn+h\r\n        do 70 it=1,10\r\n          fnew=func(xnew,yimpn,dfdy)\r\n          gy=yimp-yimpn+(9.*fnew+19.*fimp(0)-5.*fimp(1)+fimp(2))*h\/24.\r\n          dgdy=9.*h\/24.*dfdy-1.\r\n          dely=-gy\/dgdy\r\n          yimpn=yimpn+dely\r\nc         <a name=abs><font color=\"FF0000\">\r\n\t  if(abs(dely\/max(yimp,yimpn)).lt.1.e-9) go to 72\r\nc         <\/font>  ^^^^^^^^^^^^^^^^^^^^^^^^^^^\r\n   70   continue\r\n        write(6,*) 'Implicit iteration failed'\r\n        stop\r\nc\r\nc   Standard Runge-Kutta continues\r\nc\r\n   72 \txk1=h*func(xn,yrk,dfdy)\r\n      \txk2=h*func(xn+.5*h,yrk +.5*xk1,dfdy)\r\n      \txk3=h*func(xn+.5*h,yrk +.5*xk2,dfdy)\r\n      \txk4=h*func(xn+h,yrk +xk3,dfdy)\r\n        yrk =yrk +(xk1+2.*xk2+2.*xk3+xk4)\/6.\r\nc\r\nc   Do a straight Adams-Bashford integration\r\nc\r\n        yab=yam+(55.*ft(0)-59.*ft(1)+37.*ft(2)-9.*ft(3))*h\/24.\r\nc\r\nc   Do an Adams-Bashford predictor for Adams-Moulton\r\nc\r\n        yabt=yam+(55.*f(0)-59.*f(1)+37.*f(2)-9.*f(3))*h\/24.\r\n        ftnew=func(xnew,yabt,dfdy)\r\nc\r\nc   Now apply and Adams-Moulton correction step\r\nc\r\n        yam=yam+(9.*ftnew+19.*f(0)-5.*f(1)+f(2))*h\/24.\r\n        call answer(xnew,yanew)\r\n        write(10,2000)xnew,yanew,yrk,yab,yam,yimpn\r\n   80 continue\r\nc\r\n      close(10)\r\n      stop\r\n      end\r\nc\r\n      function func(x,y,dfdy)\r\n      use precision\r\n      real (rp) x,y,dfdy,func\r\nc\r\nc   A specific function f for the right hand side of the ODE\r\nc\r\n         func=1.-y**2\r\nc\r\nc   Limit the value to keep the code running when one numerical solution\r\nc   method goes nuts\r\nc\r\n         func=min(1.e10_rp,max(func,-1.e10_rp))\r\n         dfdy=-2*y\r\n      return\r\nc<a name=\"ef\"><font color=\"FF0000\">\r\n      end function func\r\nc<\/a><\/font>\r\n      subroutine answer(x,y)\r\n      use precision\r\n      real (rp) x,y,e2x\r\nc\r\nc  For the function in \"func\" this is the analytic solution to the ODE\r\nc<a name=3><font color=FF0000>\r\n      e2x=exp(2*x)\r\nc<\/font>\r\n      y=(e2x-1.)\/(1.+e2x)\r\n      return\r\n      end\r\nc<\/pre>\r\nc<\/body>\r\nc<\/html><\/pre>\n","protected":false},"excerpt":{"rendered":"
c<html> c<head><title>odeint.f<\/title><\/head> c<body> c<pre> c<a name=”module”><font color=”FF0900″> module precision c<\/font><\/a> c<a name=2><font color=FF0000> integer, parameter :: rp=selected_real_kind(13,300) c<\/font> end module precision program odeint use precision implicit none c c Demonstration of Runge-Kutta, Adams-Bashford, and Adams-Moulton for c solution of an Ordinary Differential Equation in the form: c c d y c — = f ( […]<\/p>\n","protected":false},"author":22,"featured_media":0,"parent":0,"menu_order":0,"comment_status":"closed","ping_status":"closed","template":"","meta":{"_exactmetrics_skip_tracking":false,"_exactmetrics_sitenote_active":false,"_exactmetrics_sitenote_note":"","_exactmetrics_sitenote_category":0,"footnotes":""},"categories":[],"tags":[],"class_list":["post-183","page","type-page","status-publish","hentry"],"_links":{"self":[{"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/pages\/183","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/users\/22"}],"replies":[{"embeddable":true,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/comments?post=183"}],"version-history":[{"count":1,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/pages\/183\/revisions"}],"predecessor-version":[{"id":184,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/pages\/183\/revisions\/184"}],"wp:attachment":[{"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/media?parent=183"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/categories?post=183"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.professores.uff.br\/diomarcesarlobao\/wp-json\/wp\/v2\/tags?post=183"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}