One dimensional wave equation
by Reinaldo Baretti Machin

reibaretti2004@yahoo.com

Url: www.oocities.org/reibaretti2004


c one dimensional wave equation alfa**2d^2u/dx^2= d^2u/dt^2
c alfa  = velocity // method =employed forward difference method
c Lscale  is length of string //  tscale=L/velocity
c stability imposes alfa*dt/dx<1
      dimension u(0:5000,0:5000)
      equivalence ( alfa, vel)
      real Lscale, lambda
      uexact(x,t)=sin(pi*x)*cos(2.*pi*t)
      uxcero(x)=sin(pi*x)
      dudt(x)=0.
      pi=2.*asin(1.)
      vel=2.
      Lscale=1.
      dx=.01*Lscale
      dt=.5*dx/vel
      nstep=int(Lscale/dx)
      lambda=(vel*dt/dx)
c       c final time chosen tfinal=1. ntime<10,000
      tfinal=1.
      ntime=int(tfinal/dt)
      print*,'ntime=',ntime
c boundary conditions zero amplitude at both ends x=0, x=L
      do 70 i=0,ntime
      u(0,i)=0.
      u(nstep,i)=0.
70    continue
      print*,'dt,dx,nstep,ntime=',dt,dx,nstep,ntime
c  initinial condiditons//row =1 has to be filled u1=u0+(dudt)*dt
      do 10 i=1,nstep-1
      x=dx*float(i)
      u(i,0)=uxcero(x)
      u(i,1)=uxcero(x)+dt*dudt(x)
10    continue
c the most advanced time step approx.
      do 30  l=1,ntime
      do 30  k=1,nstep-1
      u(k,l+1)=2.*u(k,l)-u(k,l-1)+lambda**2*(u(k+1,l)-2.*u(k,l)
     $ +u(k-1,l))
30    continue
       print*,'tfinal=',tfinal
       do 50 i=0,nstep,5
       x=dx*float(i)
       write(6,100)x, u(i,ntime),uexact(x,tfinal)
50    continue
100   format(2x,'x,u,uexact=',3x,e10.3,3x,e11.4,3x,e11.4)
      stop
      end
      
       ntime= 400
 dt,dx,nstep,ntime=  0.00249999994  0.00999999978 100 400
 tfinal=  1.
  x,u,uexact=    0.000E+00    0.0000E+00    0.0000E+00
  x,u,uexact=    0.500E-01    0.1564E+00    0.1564E+00
  x,u,uexact=    0.100E+00    0.3090E+00    0.3090E+00
  x,u,uexact=    0.150E+00    0.4540E+00    0.4540E+00
  x,u,uexact=    0.200E+00    0.5878E+00    0.5878E+00
  x,u,uexact=    0.250E+00    0.7071E+00    0.7071E+00
  x,u,uexact=    0.300E+00    0.8090E+00    0.8090E+00
  x,u,uexact=    0.350E+00    0.8910E+00    0.8910E+00
  x,u,uexact=    0.400E+00    0.9510E+00    0.9511E+00
  x,u,uexact=    0.450E+00    0.9877E+00    0.9877E+00
  x,u,uexact=    0.500E+00    0.1000E+01    0.1000E+01
  x,u,uexact=    0.550E+00    0.9877E+00    0.9877E+00
  x,u,uexact=    0.600E+00    0.9511E+00    0.9511E+00
  x,u,uexact=    0.650E+00    0.8910E+00    0.8910E+00
  x,u,uexact=    0.700E+00    0.8090E+00    0.8090E+00
  x,u,uexact=    0.750E+00    0.7071E+00    0.7071E+00
  x,u,uexact=    0.800E+00    0.5878E+00    0.5878E+00
  x,u,uexact=    0.850E+00    0.4540E+00    0.4540E+00
  x,u,uexact=    0.900E+00    0.3090E+00    0.3090E+00
  x,u,uexact=    0.950E+00    0.1564E+00    0.1564E+00
  x,u,uexact=    0.100E+01    0.0000E+00   -0.8742E-07