Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

  • Congratulations Chriss Miller on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

Finite Differences

Status
Not open for further replies.

Welberg123

Programmer
Joined
Jan 9, 2015
Messages
5
Location
ES
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

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:
      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
 
To compute n+1, you need n and n-1

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

Frantois Jacq
 
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
 
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:
!     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
 
Doesn't work either[sad] I'm sure it's a silly mistake that's causing the results to be incoherent, but I cannot find it
 
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
 
by the way,

whats this line:
Code:
u(i,j)=c3*u(i,j-1)+[b]alfa[/b]*(u(i+1,j-1)+u(i-1,j-1))-u(i,j-2)

alfa=alpha?
 
You should add "implicit none" at the beginning of your program to help in avoiding typos
 
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:
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?
 
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)
 
Status
Not open for further replies.

Part and Inventory Search

Sponsor

Back
Top