Usually we consider independent variables x, t. The wave equation is
, with specified initial data for 
, and
appropriate boundary conditions.
Because the coefficient matrix is constant, we can change variables,
diagonalizing, to get
where 
are the eigenvalues of the coefficient matrix, and
are combinations of 
.
We will start more simply, looking at the `one-direction' wave equation
Let us assume that a ;SPMgt; 0. We need to specify initial conditions, u(x,0)=f(x). The equation says ``u is constant along lines with slope 1/a''. That is,
A solution (ignoring boundary conditions for the moment) is
The curves dx/dt = a are called characteristic curves, and along the characteristics, du/dt=0.
How can we translate these ideas into a numerical scheme? Consider a
discritization of space into subintervals of size 
, and time increments
. Let us set notation:
. If we know 
(for example,
if n=0 we know initial data), we need to determine
for all j. We could write a difference approximation for the
time and space derivatives. It seems obvious to use
for the time derivative at point j. But what
about the x-derivative? We could use centered, forward or backward differences.
The code below tests these alternatives. Run it, and, based on the solution
form given above, give an argument for why different approximations work or fail.
program wavesolve
Implicit none
Double precision x(101),ucen(101),ubk(101),ufw(101)
Double precision dt,dx,time,pi
Integer i,Nx,nsteps
common/grid/dx,dt,nx
pi=4.d0*datan(1.d0)
Nx=101
nsteps=10
time=0.d0
dx=1.d0/(Nx-1)
dt=0.5d0*dx
do i=1,Nx
x(i)=(i-1)*dx
enddo
c INITIAL DATA define ucen, then set them all equal
time=0.d0
do i=1,nx
ucen(i)=0.d0
enddo
do i=11,31
ucen(i)=dsin(2.d0*pi*(i-11.d0)*dx/(20.d0*dx))
write(*,*)ucen(i)
enddo
do i=1,nx
ufw(i)=ucen(i)
ubk(i)=ucen(i)
enddo
open(unit=1,file='outfile')
write(1,*)"At time ",time
write(1,1001)
do i=1,nx
write(1,1011)x(i),ucen(i),ubk(i),ufw(i)
enddo
do i=1,nsteps
time=time+dt
call center(ucen)
call back(ubk)
call forward(ufw)
enddo
write(1,*)"At time ",time
write(1,1001)
do i=1,nx
write(1,1011)x(i),ucen(i),ubk(i),ufw(i)
enddo
1001 format('',3x,'X',10x,'Uc',10x,'Ub',10x,'Uf')
1011 format('',f6.4,3(2x,f10.7))
stop
end
subroutine center(u)
Implicit none
Double precision u(101),unew(101)
Double precision dt,dx
Integer i,Nx
common/grid/dx,dt,nx
do i=2,nx-1
unew(i)=u(i)-0.5d0*(dt/dx)*(u(i+1)-u(i-1))
enddo
c SET BCs -- we'll use 0 here
unew(1)=0.d0
unew(nx)=0.d0
do i=1,nx
u(i)=unew(i)
enddo
return
end
subroutine back(u)
Implicit none
Double precision u(101),unew(101)
Double precision dt,dx
Integer i,Nx
common/grid/dx,dt,nx
do i=2,nx-1
unew(i)=u(i)-(dt/dx)*(u(i)-u(i-1))
enddo
c SET BCs -- we'll use 0 here
unew(1)=0.d0
unew(nx)=0.d0
do i=1,nx
u(i)=unew(i)
enddo
return
end
subroutine forward(u)
Implicit none
Double precision u(101),unew(101)
Double precision dt,dx
Integer i,Nx
common/grid/dx,dt,nx
do i=2,nx-1
unew(i)=u(i)-(dt/dx)*(u(i+1)-u(i))
enddo
c SET BCs -- we'll use 0 here
unew(1)=0.d0
unew(nx)=0.d0
do i=1,nx
u(i)=unew(i)
enddo
return
end
Notice that in all three subroutines, I only updated interior points, and set BOTH boundary values to 0. If you look at the codes, however, the forward and backward difference codes only need 1 boundary condition. If you look at the solution u=f(x-at), it is clear that a boundary condition on the left is needed, but not on the right.
A common approach to examine numerical methods is to study their stability.
That is, consider a numerical approximation
to be composed of an exact solution plus an error. If we were to advance the
approximate solution by one timestep using the algorithm, what would happen
to the error term? Stability means the error does not grow. For linear
equations, we can use Fourier methods to investigate stability.
Define the Fourier Transform
and the inverse Transform
Now, for instance, lets look at the forward time, centered space
differencing introduced in the code above. A finite difference like
introduces a 
on the Fourier side
so we get
By trig manipulations, with 
,
The multiplicative factor on the right is called the amplification factor. It describes how the numerical approximation (and any error!) evolves in time.
We define a few concepts.
A finite difference scheme is convergent if
, as 
.
A finite difference scheme 
is consistent
with a PDE L u if for a smooth
function u(x,t), 
as
.
A finite difference scheme is stable if there exists numbers
C, 
such that, for a fixed positive time T,
for 
Lax-Wendroff Theorem A consistent finite difference scheme is convergent if and only if it is stable.
Theorem A finite difference scheme (for a homogeneous PDE)
is stable iff for 
small enough,
.
For the FTCS scheme, 
,
and, taking modulus, 
.
Thus the scheme is unconditionally unstable.
For the backwards differencing formula,
, and its modulus
is less than 1 if 
.
The Lax-Friedrichs scheme overcomes the FTSC instability by replacing
on the right by its nearest-neighbor average
The scheme below is a second order modification of the Lax-Friedrichs scheme,
for a system of equations
, with boundary conditions that reflect the wave.
Homework As it appears below, the scheme runs until a time t=1.0. At that time, the exact solution ought to equal the initial data. See how well it does. Change the final time - making it longer (say, tfinal=2, to test how well the square data is preserved, and shorter (maybe tfinal=0.75), to see the waves reflecting off the boundary. Change the initial data to be a sine wave with the same initial support. How well is the sine data preserved after t=1? Change the frequency of that sine data, and tell how the preservation of the data corrolates with the frequency.
C ***********************************************************
program wavesolve
implicit double precision(a-h,o-z)
c SOLVE LINEAR WAVE SYSTEM U_t + A U_x = 0
c U=(u,v), A=[ [0 -1], [-1 0] ]
c USING TADMOR's 2nd-ORDER LAX-FRIEDRICHS SCHEME
c
dimension U(2,102),x(102),Ux(2,102),Uhalf(2,102)
common/grid/dt,dx,dtdx,nx
tfinal=1.0d0
nx=102
dx=1.d0/(nx-2)
do i=1,nx
x(i)=(i-1-0.5)*dx
enddo
c INITIAL DATA
time=0.d0
do i=1,nx
do j=1,2
u(j,i)=0.d0
enddo
enddo
do i=4*nx/10+2,6*nx/10
u(1,i)=1.d0
enddo
call print(U,x,time)
c MAIN ITERATION LOOP
dt=0.9*dx
iend=0
do iter=1,100
do inner=1,5
tnew=time+dt
if (tnew .gt. tfinal) then
dt=tfinal-time
tnew=time+dt
iend=1
endif
dtdx=dt/dx
call slope(U,Ux)
call midtime(U,Ux,Uhalf)
call update(U,Ux,Uhalf)
time=tnew
if (iend .eq. 1) then
call print(U,x,time)
go to 999
endif
enddo
call print(U,x,time)
enddo
999 stop
end
C ***********************************************************
C ***********************************************************
subroutine slope(U,Ux)
implicit double precision(a-h,o-z)
dimension U(2,102),Ux(2,102)
common/grid/dt,dx,dtdx,nx
c STATEMENT FUNCTION sign
sgn(w)=1.d0*DSIGN(1.d0,w)
gamma=2.d0
do i=1,nx
Ux(1,i)=0.d0
Ux(2,i)=0.d0
enddo
do i=2,nx-1
df=U(1,i+1)-U(1,i)
db=U(1,i)-U(1,i-1)
df=df+db
Ux(1,i)=0.5d0*(sgn(df)+sgn(db))*
1 dmin1(gamma*dabs(df),gamma*dabs(db),0.5d0*dabs(dc))
df=U(2,i+1)-U(2,i)
db=U(2,i)-U(2,i-1)
df=df+db
Ux(2,i)=0.5d0*(sgn(df)+sgn(db))*
1 dmin1(gamma*dabs(df),gamma*dabs(db),0.5d0*dabs(dc))
enddo
return
end
C ***********************************************************
C ***********************************************************
subroutine midtime(U,Ux,Uhalf)
implicit double precision(a-h,o-z)
dimension U(2,102),Ux(2,102),Uhalf(2,102)
common/grid/dt,dx,dtdx,nx
do i=2,nx-1
Uhalf(1,i)=U(1,i) +0.5d0*dtdx*Ux(2,i)
Uhalf(2,i)=U(2,i) +0.5d0*dtdx*Ux(1,i)
enddo
call BCs(Uhalf)
return
end
C ***********************************************************
C ***********************************************************
subroutine update(U,Ux,Uhalf)
implicit double precision(a-h,o-z)
dimension U(2,102),Ux(2,102),Uhalf(2,102),Unew(2,102)
common/grid/dt,dx,dtdx,nx
do i=2,nx-1
Unew(1,i)=0.5d0*(U(1,i+1)+U(1,i-1))
1 -0.25d0*(Ux(1,i+1)-Ux(1,i-1))
2 +0.5d0*dtdx*(Uhalf(2,i+1)-Uhalf(2,i-1))
Unew(2,i)=0.5d0*(U(2,i+1)+U(2,i-1))
1 -0.25d0*(Ux(2,i+1)-Ux(2,i-1))
2 +0.5d0*dtdx*(Uhalf(1,i+1)-Uhalf(1,i-1))
enddo
call BCs(Unew)
do i=1,nx
do j=1,2
U(j,i)=Unew(j,i)
enddo
enddo
return
end
C ***********************************************************
C ***********************************************************
subroutine BCs(Us)
implicit double precision(a-h,o-z)
dimension Us(2,102)
common/grid/dt,dx,dtdx,nx
C u reflection, v inverse reflection
Us(1,1)=Us(1,2)
Us(1,nx)=Us(1,nx-1)
Us(2,1)=-Us(2,2)
Us(2,nx)=-Us(2,nx-1)
return
end
C ***********************************************************
C ***********************************************************
subroutine print(U,x,time)
implicit double precision(a-h,o-z)
dimension U(2,102),x(102)
common/grid/dt,dx,dtdx,nx
write(*,1000)time
write(*,1001)
do i=1,nx
write(*,1010)x(i),U(1,i),U(2,i)
enddo
1000 format('','AT TIME ',f9.7)
1001 format('',3x,'X',12x,'u',12x,'v')
1010 format('',f9.6,2(2x,f10.7))
return
end
C ***********************************************************