How to Modify a Finite Difference Method Code for a Non-Zero Boundary Condition?

  • Thread starter Thread starter fahad20
  • Start date Start date
fahad20
Messages
7
Reaction score
0
hello...i m trying to solve this ordinary differential equation problem..but it is giving errors..pleasezz help

problem: d^2 T/Dr^2 + (1/r) DT/Dr=1+ (r/R)^2

T'(0)=0 and T(R)=0


solve this problem for T(r),where R=1.0

SOLUTION I HAVE CUM UP WITH:

put R=1 in this equation...

now let v=T'
T''=V'

so

DT/Dr=v

Dv/dr=1+(T**2/V**2)-(V**2/T)

V(0)=0 and T(R)=0

and now here's the code...



dimension p(100,4),r(100),xn(100),xt(100)

write(*,10)

10 format(1x,'input number of equation:')

read(*,*)N
write(*,12)

12 format(1x,'enter the initial value of the independent variable:')

read(*,*)T

write(*,14)

14 format(1x,'enter the increment of the independent variable:')

read(*,*)dt

write(*,16)

16 format(1x,'enter the ending value of the independent variable:')

read(*,*)tend

write(*,18)

18 format(1x,'enter the initial value of the dependent variables:')

do 20 i=1,n

20 read(*,*)r(i)

write(*,22)

22 format(9x,'t',5x,'r(1) through r(n)')

write(*,24) T,(r(i),i=1,n)

24 format(f10.3,5e14.5)

26 call rkn(x,n,dt,t,p,xt,xn)

t=t+dt

do 30 i=1,n

30 r(i)=xn(i)

write(*,24)t,(r(i),i=1,n)

if (t.ge.tend) goto 50

goto 26

50 write(*,52)

52 format(1x,'r-k parameters are:')

do 54 i=1,n

54 write(*,56) (p(i,j),j=1,4)

56 format(4e15.5)

end

function f(r,t,i,n)
dimension r(n)

goto (10,20),i

10 f=r(1)

return

20 f=1+(r(2)**2/r(1)**2)-(r(1)**2/r(2))


return

end


subroutine rkn(xin,n,dt,t,p,xt,xout)

! solving n first-order ode d(xin)/dt=f(xin,t)

!knowing xin at t,this sub.finds xout at t+dt

!p are the runge-kutta 4th-order parameters
dimension p(n,4),xin(n),xout(n),xt(n)

do 5 i=1,n

5 p(i,1)=dt*f(xin,t,i,n)

do 15 j=2,3

do 10 i=1,n

10 xt(i)=xin(i)+p(i,j-1)/2.

do 15 i=1,n

15 p(i,j)=dt*f(xt,t+dt/2.,i,n)

do 20 i=1,n

20 xt(i)=xin(i)+p(i,3)

do 25 i=1,n

p(i,4)=dt*f(xt,t+dt,i,n)


25 xout(i)=xin(i)+(p(i,1)+2.*(p(i,2)+p(i,3))+p(i,4))/6.

return

end




but it is giving eror..please correct this..thank u
 
Physics news on Phys.org


You made several errors, I have bolded them for your below.

fahad20 said:
hello...i m[/b trying to solve this ordinary differential equation problem..but it is giving errors..pleasezz help

problem: d^2 T/Dr^2 + (1/r) DT/Dr=1+ (r/R)^2

T'(0)=0 and T(R)=0solve this problem for T(r),where R=1.0

SOLUTION I HAVE CUM UP WITH:

put R=1 in this equation...

now let v=T'
T''=V'

so

DT/Dr=v

Dv/dr=1+(T**2/V**2)-(V**2/T)

V(0)=0 and T(R)=0

and now heres the code...
...but it is giving eror..please correct this..thank u


I didn't bother checking sentence structure or punctuation. If you'd like, I can correct those for you as well.
 


fahad20 said:
but it is giving errors..

What kind of errors?

If the program doesn't even compile, what error messages does the compiler give you?

If the program compiles but doesn't run correctly, show us

(a) some sample input

(b) the output that the program should produce

(c) the output that the program actually produces
 


fahad20 said:

Homework Statement



solve for T(r)

Homework Equations





The Attempt at a Solution



hello...i m trying to solve this ordinary differential equation problem..but it is giving errors..pleasezz help

problem: d^2 T/Dr^2 + (1/r) DT/Dr=1+ (r/R)^2

T'(0)=0 and T(R)=0


solve this problem for T(r),where R=1.0

SOLUTION I HAVE CUM UP WITH:

put R=1 in this equation...

now let v=T'
T''=V'

so

DT/Dr=v

Dv/dr=1+(T**2/V**2)-(V**2/T)

V(0)=0 and T(R)=0

and now here's the code...



dimension p(100,4),r(100),xn(100),xt(100)

write(*,10)

10 format(1x,'input number of equation:')

read(*,*)N
write(*,12)

12 format(1x,'enter the initial value of the independent variable:')

read(*,*)T

write(*,14)

14 format(1x,'enter the increment of the independent variable:')

read(*,*)dt

write(*,16)

16 format(1x,'enter the ending value of the independent variable:')

read(*,*)tend

write(*,18)

18 format(1x,'enter the initial value of the dependent variables:')

do 20 i=1,n

20 read(*,*)r(i)

write(*,22)

22 format(9x,'t',5x,'r(1) through r(n)')

write(*,24) T,(r(i),i=1,n)

24 format(f10.3,5e14.5)

26 call rkn(x,n,dt,t,p,xt,xn)

t=t+dt

do 30 i=1,n

30 r(i)=xn(i)

write(*,24)t,(r(i),i=1,n)

if (t.ge.tend) goto 50

goto 26

50 write(*,52)

52 format(1x,'r-k parameters are:')

do 54 i=1,n

54 write(*,56) (p(i,j),j=1,4)

56 format(4e15.5)

end

function f(r,t,i,n)
dimension r(n)

goto (10,20),i

10 f=r(1)

return

20 f=1+(r(2)**2/r(1)**2)-(r(1)**2/r(2))


return

end


subroutine rkn(xin,n,dt,t,p,xt,xout)

! solving n first-order ode d(xin)/dt=f(xin,t)

!knowing xin at t,this sub.finds xout at t+dt

!p are the runge-kutta 4th-order parameters
dimension p(n,4),xin(n),xout(n),xt(n)

do 5 i=1,n

5 p(i,1)=dt*f(xin,t,i,n)

do 15 j=2,3

do 10 i=1,n

10 xt(i)=xin(i)+p(i,j-1)/2.

do 15 i=1,n

15 p(i,j)=dt*f(xt,t+dt/2.,i,n)

do 20 i=1,n

20 xt(i)=xin(i)+p(i,3)

do 25 i=1,n

p(i,4)=dt*f(xt,t+dt,i,n)


25 xout(i)=xin(i)+(p(i,1)+2.*(p(i,2)+p(i,3))+p(i,4))/6.

return

end




but it is giving eror..please correct this..thank u

You have several do loops that don't have continue statements. The loops should look something like this:
Code: 
do 10 i = 1, n
     sum = sum + i
     write(*,*) 'i =', i
     write(*,*) 'sum =', sum
10  continue

Your code is very hard to read with statement labels on input and output statements that you probably don't need, and with gotos here and there.
 


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...





!Fortran Version

! Program OdeBcpFD – Ordinary differential equation, boundry-value Problem

! 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...

thank u ..
 


Is this a programming course? Are you required to use that program?

If, after taking v= T', you multiply the original equation, v'+ v/r= 1+ r^2 by r, you get rv'+ v= (rv)'= r+ r^3 and can solve for v just by integrating both sides:
rv= (1/2)r^2+ (1/4)r^4+ C so that v= T'= (1/2)r+ (1/4)r^3+ C/r and T can be found by integrating again.
 


please tryy to understnd...

the code i have is for ..

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

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

(both B.C's are ZERO)

this is a finite difference method code in fortran..


now i have this problem...

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

in which the B.C's are not zero... ----> T(1)100 and T(2)=0

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

how will this code change for above problem..
 

Similar threads

Direct Proof that every zero of p(T) is an eigenvalue of T
Replies
24
Views
3K
Avoid unpleasant integrals in solving IVP
Replies
3
Views
2K
Integrating vector-valued functions along curves
Replies
5
Views
1K
Dynamical Systems: how to find equation for Poincare map?
Replies
1
Views
2K
Use Stokes' theorem on intersection of two surfaces
Replies
3
Views
2K
How Do Vector Spaces of Linear Maps Differ from Standard Vector Spaces?
Replies
9
Views
2K
Estimator Exercise Homework Solution
Replies
11
Views
2K
Is the Calculation of the Vector Line Integral Over a Square Correct?
Replies
12
Views
2K
Mean and var of an exponential distribution using Fourier transforms
Replies
2
Views
2K
A Show positivity and boundedness of a non-linear system
Replies
2
Views
484

Hot Threads

Recent Insights

Back
Top