Coupled Differential Equations


Nusc
06-09-2009, 11:54 PM
Hi,

does anyone here know where I can find an algorithm to solve coupled differential equations with FORTRAN ?

Thanks,

wazza129
12-27-2010, 06:55 AM
a lill help plzzz.....i have this code for a 2nd ODE..but with boundary conditions.

this code is for ..

d^2 Z/dR^2 + (1/Z) dZ/dR=-p/T

Z(Ri)=0 and Z(Ro)=0

this is a finite difference method code in fortran..

please c...







!

! Solved by finite-difference method . UP to 99 stations.

DIMENSION X(99),C(99,99),R(99)

WRITE (*,2)

2 FORMAT(1X,' Program OdeBvpFD finite difference solution of O.D.E. Boundary-Value Problem.'/)

Write (*,4)

4 Format(1X,'have you Edited The functions CIJ And RI for generating the elements of [C] and {R} in [C] {Y}={R}?')

READ (*,*) I1

Write (*,6)

6 format(1x,'enter the number of (in-between,excluding boundaries)stations and stepsize:')


read(*,*)n,dx

write(*,8)
8 Format (1X,'Enter the First (left Boundary) X Value :' )

READ (*,*) X1

! Calculate [C] and {R}

Do 10 I=1,N

X (I) =X1+I*DX

R (I) =RI(x(i),dx)

DO 10 J=1,N

10 c(i,j)=cij(x(i),i,j,dx)

CALL GAUSS (C,N,99,R)




! {Y} is in {R}.

WRITE (*,12) (R(k),K=1,N)

12 Format(1x,'The Solution is : '//5(5E16.5/))

STOP

END

Function RI (X,DX)

RI=-5*DX **2 /100.

Return

End

FUNCTION CIJ (X,I,J,DX)

CIJ=100.

IF (I.EQ.j) CIJ=-2.

IF (I.EQ.(j-1)) CIJ=1.-DX/2./X

IF (I.EQ.(J+1) ) CIJ = 1.+DX/2./X

RETURN

END

SUBROUTINE GAUSS ( C , N, M, R)



!

!SOLVES MATRIX EQUATION C(N,N) *X(N)=R(N) BY GAUSSIAN ELIMINATION.

!X and V share same storage. C is dimensioned (M,M) in the calling program.

!

DIMENSION C(M,M),r(N)

!

N1=N-1

DO 25 K=1,N1

KP1=K+1

!

!NORMALIZATION

!

DO 10 J=KP1,N

10 C(K,J)=C(K,J)/C(K,K)

r(K)=r(K)/C(K,K)

!

!ELIMINATION

!

DO 25 I=KP1,N

DO 20 J=KP1,N

20 C(I,J)=C(I,J)-C(I,K)*C(K,J)

25 r(I)=r(I)-C(I,K)*r(K)

!

!BACKWARD SUBSTITUTION

!

r(N)=r(N)/C(N,N)

DO 35 I=1 ,N1

IR=N-I

IR1=IR+1

DO 35 J=IR1,N

35 r(IR)=r(IR)-C(IR,J)*r(J)

RETURN

END



now i have this problem...

d^2 T/dr^2 + (1/r) dT/dr=0

and B.C's are T(1)100 and T(2)=0

how will this code change..especially ..how to cater for these two B.Cs with r non-zero...

i.e. at r=1,T=100 and at r=2,T=0

thank u ..

a prompt reply needed