Finite Differences

[Welberg123]

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:
d[sup]2[/sup]V(x,t)/dt[sup]2[/sup] + 2b·dV(x,t)/dt - v[sup]2[/sup]·d[sup]2[/sup]V(x,t)/dx[sup]2[/sup] = 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[sup](-bt)[/sup], and the variables τ=bt and z=bx/v, then I get:

d[sup]2[/sup]u(z,τ)/dτ[sup]2[/sup] - u(z,τ) - d[sup]2[/sup]u(z,τ)/dz[sup]2[/sup]= 0

And translating into Finite Differences: u[sup]n+1[/sup][sub]j[/sub]= [2+(Δτ)[sup]2[/sup]-2α]u[sup]n[/sup][sub]j[/sub]+α(u[sup]n[/sup][sub]j+1[/sub]+u[sup]n[/sup][sub]j-1[/sub])-u[sup]n-1[/sup][sub]j[/sub]

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

Code:

Spoiler

      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

 

[FJacq]

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

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

Frantois Jacq

 

[Welberg123]

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

 

[FJacq]

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

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

You could try :

!     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

 

[Welberg123]

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

 

[Welberg123]

Doesn't work either

 

[pingat]

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

 

[pingat]

by the way,

whats this line:

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?

 

[salgerman]

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

 

[Welberg123]

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:

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?

 

[salgerman]

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)

 

[salgerman]

Sorry...bali-->valid