!******************************************************************** !* Differential equations with p variables of order 1 * !* by Runge-Kutta method of order 4 * !* ---------------------------------------------------------------- * !* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 by Marc * !* DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991" * !* * !* F90 version 1.1 By J-P Moreau, Paris * !* (www.jpmoreau.fr) * !* ---------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* Example #1: integrate system of equations from x=0 to x=3: * !* y1' = y2 + y3 - 3*y1 * !* y2' = y1 + y3 - 3*y2 * !* y3' = y1 + y2 - 3*y3 * !* with initial conditions: y1(0)=1, y2(0)=2 and y3(0)=-1 * !* * !* * !* DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1 * !* of type yi' = f(y1,y2,...,yn), i=1..n * !* * !* number of variables: 3 * !* begin value x : 0 * !* end value x : 3 * !* y1 value at x0 : 1 * !* y2 value at x0 : 2 * !* y3 value at x0 : -1 * !* number of points : 7 * !* finesse : 30 * !* * !* X Y1 Y2 Y3 * !* -------------------------------------------------------- * !* 0.000000 1.000000 2.000000 -1.000000 * !* 0.500000 0.449466 0.584801 0.178795 * !* 1.000000 0.251358 0.269674 0.214727 * !* 1.500000 0.149580 0.152058 0.144622 * !* 2.000000 0.090335 0.090671 0.089664 * !* 2.500000 0.054738 0.054784 0.054648 * !* 3.000000 0.033193 0.033200 0.033181 * !* -------------------------------------------------------- * !* * !* Example #2: integrate system of equations from x=0 to PI*Sqrt(2) * !* y1' = y2 * !* y2' = -4*y1 - 3*y3 * !* y3' = y4 * !* y4' = -8*y1 - 2*y3 * !* with initial conditions: y1(0)=3, y2(0)=0, y3(0)=4, y4(0)=0 * !* * !* * !* DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1 * !* of type yi' = f(y1,y2,...,yn), i=1..n * !* * !* number of variables: 4 * !* begin value x : 0 * !* end value x : 4.442883 * !* y1 value at x0 : 3 * !* y2 value at x0 : 0 * !* y3 value at x0 : 4 * !* y4 value at x0 : 0 * !* number of points : 9 * !* finesse : 30 * !* * !* X Y1 Y2 Y3 Y4 * !* -------------------------------------------------------- * !* 0.000000 3.000000 0.000000 4.000000 0.000000 * !* 0.555360 0.000000 -8.485281 0.000000 -11.313708 * !* 1.110721 -3.000000 -0.000001 -4.000000 -0.000002 * !* 1.666081 -0.000001 -8.485281 -0.000001 11.313708 * !* 2.221442 3.000000 0.000003 4.000000 0.000003 * !* 2.776802 0.000001 -8.485281 0.000002 -11.313708 * !* 3.332162 -3.000000 -0.000004 -4.000000 -0.000005 * !* 3.887523 -0.000002 8.485281 -0.000002 11.313708 * !* 4.442883 3.000000 0.000005 4.000000 0.000007 * !* -------------------------------------------------------- * !* * !* Release 1.1: Added example #2 (Sept. 2008). * !******************************************************************** PROGRAM TEST_EQDIFP integer fi,i,ndata,p real*8 xi,xf real*8 yi(10) print *,' ' print *,' DIFFERENTIAL EQUATIONS WITH P VARIABLES OF ORDER 1' print *,' of type yi'' = f(y1,y2,...,yn), i=1..n' print *,' ' write(*,"(' number of variables : ')",advance='no') read *, p print *,' ' write(*,"(' begin value x : ')",advance='no') read *, xi write(*,"(' end value x : ')",advance='no') read *, xf do i=0, p-1 write(*,50,advance='no') i+1 read *, yi(i) end do write(*,"(' number of points : ')",advance='no') read *, ndata write(*,"(' finesse : ')",advance='no') read *, fi call Eqdifp(xi,xf,yi,ndata,p,fi); 50 format(' y',i1,' value at x0 : ') end !Example #1: y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3 !real*8 function fp(k, x, y) ! integer k; real*8 x, y(10) ! if (k.eq.0) fp = y(1)+y(2)-3.d0*y(0) ! if (k.eq.1) fp = y(0)+y(2)-3.d0*y(1) ! if (k.eq.2) fp = y(0)+y(1)-3.d0*y(2) !end !Exemple #2: y1'=y2, y2'=-4y1-3y3, y'3=y4, y'4=-8y1-2y3 } real*8 function fp(k, x, y) integer k real*8 x, y(10) if (k==0) fp=y(1) if (k==1) fp=-4.d0*y(0)-3.d0*y(2) if (k==2) fp=y(3) if (k==3) fp=-8.d0*y(0)-2.d0*y(2) End ! print results Subroutine Display(t1,t2,t3,t4,m,p,xi,xf) parameter(SIZE=50) real*8 t1(SIZE),t2(SIZE),t3(SIZE),t4(SIZE),xi,xf integer i,m,p real*8 h,x h=(xf-xi)/(m-1) x=xi-h print *,' ' write(*,90,advance='no') do i=1, p write(*,91,advance='no') i end do print *,' ' print *,'------------------------------------------------------' do i=1, m x=x+h if (p==2) then write(*,100) x,t1(i),t2(i) else if (p==3) then write(*,101) x,t1(i),t2(i),t3(i) else if (p > 3) then write(*,102) x,t1(i),t2(i),t3(i),t4(i) end if end do print *,'------------------------------------------------------' print *,' ' return 90 format(' X') 91 format(' Y',I1) 100 format(' ',F9.6,' ',F10.6,' ',F10.6) 101 format(' ',F9.6,' ',F10.6,' ',F10.6,' ',F10.6) 102 format(' ',F9.6,' ',F10.6,' ',F10.6,' ',F10.6,' ',F10.6) end !**************************************************************************** !* SOLVING DIFFERENTIAL SYSTEMS WITH P VARIABLES OF ORDER 1 * !* of type yi' = f(y1,y2,...,yn), i=1..n * !* ------------------------------------------------------------------------ * !* INPUTS: * !* xi, xf begin, end values of variable x * !* Yi table of begin values of functions at xi * !* m number of points to display * !* p number of independant variables * !* fi finesse (number of intermediary points) * !* * !* OUTPUTS: * !* t1,t2 real vectors storing the results for first two functions, * !* y1 and y2 (passed to subroutine Affiche1(). * !* ------------------------------------------------------------------------ * !* EXAMPLE: y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3 * !* Exact solution : y1 = 1/3 (exp(-4x) + 2 exp(-x)) * !* y2 = 1/3 (4exp(-4x) + 2 exp(-x)) * !* y3 = 1/3 (-5exp(-4x)+ 2 exp(-x)) * !**************************************************************************** Subroutine Eqdifp(xi,xf,Yi,m,p,fi) parameter(SIZE=50) real*8 t1(SIZE),t2(SIZE),t3(SIZE),t4(SIZE),xi,xf,Yi(10) integer fi,m,p real*8 h,x,fp real*8 ta(10),tb(10),tc(10),td(10),y(10),z(10) integer i,j,k,ni if (fi<1) return h = (xf - xi) / fi / (m-1) p=p-1 t1(1)=Yi(0) t2(1)=Yi(1) t3(1)=Yi(2) t4(1)=Yi(3) do k=0, p y(k)=Yi(k); z(k)=Yi(k) end do do i=1, m ni=(i-1)*fi-1 do j=1, fi x=xi+h*(ni+j) do k=0, p y(k)=z(k) end do do k=0, p ta(k)=h*fp(k,x,y) end do do k=0, p y(k)=z(k)+ta(k)/2.d0 end do x=x+h/2.d0 do k=0, p tb(k)=h*fp(k,x,y) end do do k=0, p y(k)=z(k)+tb(k)/2.d0 end do do k=0, p tc(k)=h*fp(k,x,y) end do do k=0, p y(k)=z(k)+tc(k) end do x=x+h/2.d0 do k=0, p td(k)=h*fp(k,x,y) end do do k=0, p z(k)=z(k)+(ta(k)+2.d0*tb(k)+2.d0*tc(k)+td(k))/6.d0 end do end do t1(i+1)=z(0) t2(i+1)=z(1) t3(i+1)=z(2) t4(i+1)=z(3) end do call Display(t1,t2,t3,t4,m,p+1,xi,xf) return end ! End of file teqdifp.f90