INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

#### Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

# Finite Differences

## Finite Differences

(OP)
Hi,

I'm trying to write a program in Fortran77 to solve the telegraph's equation using Finite Differences, but there seems to be a problem somewhere in the code. This is what I got so far:

Explanation of the idea

#### Spoiler:

Original equation:
d2V(x,t)/dt2 + 2b·dV(x,t)/dt - v2·d2V(x,t)/dx2 = 0

where b-1 is a time parameter and v the propagation speed. Now, to make the equation user-friendly, I used V(x,t)=u(x,t)·e(-bt), and the variables τ=bt and z=bx/v, then I get:

d2u(z,τ)/dτ2 - u(z,τ) - d2u(z,τ)/dz2= 0

And translating into Finite Differences: un+1j= [2+(Δτ)2-2α]unj+α(unj+1+unj-1)-un-1j

where α=((Δτ)2)/((Δz)2). The idea is to solve this equation with the Initial condition V(z,τ=0)=C1·e(-C2·z2)

Code:

#### CODE --> fortran77

program Telegraph

integer*4 Nspatial, Ntemporal, Npoints, i,j
real*8 dx, dt, p, L, T(1001,10000), u(1001,10000), c1, alpha,c3,c2

open(unit=1, file="Telegraph.txt")

!     Variables:

Ntemporal= 500
dx=0.001
dt = 0.001
Npoints=100
alpha=(dt*dt)/(dx*dx)
C1=0.7283656
C2=5.5555555
C3=(2.0d0+2.0d0*dt-2.0d0*alpha)

!!!!! Main program:

write(1,*) "---------------------------------------------------"
write(1,*) "    u(i,j)   ||     time      ||    position       "
write(1,*) "---------------------------------------------------"

!     Initial conditions:
do i=1, Npoints-1
u(i,1)=C1*exp(-C2*(i**2.0d0))
enddo

!     Iteration:
do j=3, Ntemporal
do i=2, Npoints
u(i,j)=c3*u(i,j-1)+alfa*(u(i+1,j-1)+u(i-1,j-1))-u(i,j-2)
write(1,*) u(i,j), i, j
end do
end do

write(*,*) "----------------------------------------"
write(*,*) "u(i,j) saved at Telegraph.txt"
write(*,*) "----------------------------------------"

end 

As you can see, the result is nonsensical (I guess due to some mistake in the iteration), if someone can pinpoint where have I gone wrong, it would be of great help.

Thanks

### RE: Finite Differences

To compute n+1, you need n and n-1

=> you need to initialize u for n=1 AND for n=2 !

Frantois Jacq

### RE: Finite Differences

(OP)
Thanks for the answer. If I start the temporal loop (j index) at j=1, the result still doesn't make sense. I also tried to add a reflexion boundary condition to avoid those two values:
do j=1, 2
u(i,j)=u(i,-j)
enddo

but still doesn't work. Could you please be more specific? Thanks again

### RE: Finite Differences

do j=1, 2
u(i,j)=u(i,-j)
enddo

But you cannot write that ! What means u(i,-j) ????

You could try :

#### CODE --> fortran

!     Initial conditions:
do i=1, Npoints-1
u(i,1)=C1*exp(-C2*(i**2.0d0))
u(i,2)=u(i,1)
enddo 

Frantois Jacq

### RE: Finite Differences

(OP)
Doesn't work either I'm sure it's a silly mistake that's causing the results to be incoherent, but I cannot find it

### RE: Finite Differences

(OP)
Doesn't work either

### RE: Finite Differences

Hey there :)

I guess its a question of stability. Haven't looked into the physics, but dx=dt may causes problems.
try a bigger dx and a smaller dt.

greetings from austria

### RE: Finite Differences

by the way,

whats this line:

#### CODE --> fortran

u(i,j)=c3*u(i,j-1)+alfa*(u(i+1,j-1)+u(i-1,j-1))-u(i,j-2)

alfa=alpha?

### RE: Finite Differences

You should add "implicit none" at the beginning of your program to help in avoiding typos

### RE: Finite Differences

(OP)
Thanks all for the answers, *alfa was a typo when I wrote it here, in my program it was written correctly. I've been looking deeper into the physics and I made a couple of mistakes with the boundary conditions, here's the updated code:

#### CODE --> fortran77

program Telegraph

integer*4 Nspatial, Ntemporal, Npoints, i,j
real*8 dx, dt, p, L, T(1001,10000), u(1001,10000), c1, alpha,c3,c2

open(unit=1, file="Telegraph.txt")
open(unit=2, file="u(i,j).txt")

!     Variables:

Ntemporal= 500
dx=0.01
dt=0.001
Npoints=20
alpha=(dt*dt)/(dx*dx)
C1=0.7283656
C2=5.5555555
C3=(2.0d0+2.0d0*dt-2.0d0*alpha)

write(1,*) "---------------------------------------------------"
write(1,*) "    u(i,j)   ||     time      ||    position       "
write(1,*) "---------------------------------------------------"

!!!!! Initial conditions:
do i=1, Npoints
u(i,0)=C1*exp(-C2*(i**2.0d0))
enddo

do i=1, Npoints
u(i,1)=0.50d0*u(i,0)*(C3+2.0d0*dt)+(alpha/2)*(u(i+1,0)+u(i-1,0))
enddo

!!!!! Loop:
do j=2, Ntemporal
do i=1, Npoints
u(i,j)=C3*u(i,j-1)+alpha*(u(i+1,j-1)+u(i-1,j-1))-u(i,j-2)
end do
end do

do j=2, Ntemporal
do i=1,Npunts
x=dx*dble(i-1)
time=dt*dble(j-1)
write(1,100) x,time
write(2,*) u(i,j)
100   format(F12.3,F15.3)
enddo
enddo

write(*,*) "----------------------------------------"
write(*,*) "u(i,j) saved at Telegraph.txt"
write(*,*) "----------------------------------------"

end 

I'm pretty sure the set up now is the right one, but I get the famous warning message: "Warning: Array reference at (1) is out of bounds (0 < 1) in dimension 2", where (1) is u(i,0). Any ideas as how to solve this?

### RE: Finite Differences

Zero is not a Bali index value... Specifying a single value defines 1as the default first index... If you want to start at zero, you need to specify the range explicitly... U(1000,0:10000)

### RE: Finite Differences

Sorry...bali-->valid

#### Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

#### Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

Close Box

# Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!