INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Log In

Come Join Us!

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • 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.

Jobs

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:

Spoiler:

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

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.

Reply To This Thread

Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members!

Resources

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:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close