# Implementing Explicit formulation of 1D wave equation in Matlab

So the theory is straightforward. We have: $$\frac{\partial^2U}{\partial t^2}=c^2 \frac{\partial^2U}{\partial x^2}$$

discretizing it gives:

$$\frac{U(i+1,j)- 2U(i,j) + U(i-1,j)}{(\Delta t)^2} = c^2 \frac{U(i,j+1)-2U(i,j) + U(i,j-1)}{(\Delta x)^2}$$

where $U(i,j)$ is $U(t_i,x_j)$ with $t_i=i*\Delta t$ and $x_j=j*\Delta x$;

Now, coding it up in Matlab I have:

    function [u_new] = Explicit_W(dx,dt,n,initialize_pos,initialize_vel)
%Inputs
X=1;   %X: length of the space domain
J=X/dx; %J: Total # of length intervals
f=str2func(initialize_pos); %initial conditions: u(0,x)
g=str2func(initialize_vel); %initial conditions: du/dt(0,x)
s2=(dt/dx)^  2;
c=2;
%Initializing the system,
%u_ancient, u_old and u_new are U(i-1,:),U(i-1,:) and U(i+1,:) respectively.
u_ancient=f(J);
u_old=    f(J)+ dt*g(J)+ 1/2*c^2*Tridiag(J+1,s2,-2*s2,s2)*f(J);
u_new=    Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_old  - u_ancient;

for i=1:n-2, % iterate n times, 2 of which are done already!
%Iterating the system using the Explicit formulation
temp=u_new;
u_new=Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_old - u_ancient;
u_ancient=u_old;
u_old=temp;
end

end


where Tridiag returns a matrix with entries on the tri-diagonals and its defined to be:

function [ A ] = Tridiag( N,r1,r2,r3 )
%This function creates a sparse tridiagonal matrix with fixed values
%on the Diagonals.
%   comes in handy when solving PDEs!
i=[1:N 2:N   1:N-1]; %row indices
j=[1:N 1:N-1 2:N  ]; %column indices
s=[(r2)*ones(1,N) r1*ones(1,N-1) r3*ones(1,N-1)]; %nonzero elements
A=sparse(i,j,s); % create the matrix
end


and for f, g I use Sin2 and Zero respectively, which are defined to be:

function [ u0 ] = Sin2(L)
%TENT Summary of this function goes here
%   Initializes the array to Sin(2*pi*x)
u0=zeros(floor(L+1),1);
X=1;
x=(0:X/L:1)';
for j=1:(L+1),
u0=sin(2*pi*x);
end
end


and

function [ u0 ] = Zero( L)
%Summary of this function goes here
%   Initializes the array to a zero!
u0=zeros(floor(L+1),1);

end


Now, solving it analytically, if use f=sin2 and g=zeros, we know the answer should be:

$$u(x,t)=\frac{\sin(2\pi(x+ct))+\sin(2\pi(x-ct)}{2}$$

so the expected range is $\{-1,1\}$. but when I execute Explicit_W, say:

Explicit_W(0.05,0.0010,30,'Sin2','Zero')


I get funny results, ranging from -124 to 124! Here is a picture of the results: http://i.stack.imgur.com/wON52.png

So it has the form of the solution, but something is scaling it up! I have checked the sub-functions used in this on first order problems so I know there is nothing wrong with them. Whatever is causing the problem is in the main function. Any idea why it's acting the way it does?

UPDATE:

So it was brought to my attention that I had forgotten to implement the boundary conditions, here $u(0,t)=u(1,t)=0.$ Silly mistake! So the modified section of the code is:

u_ancient=f(J);
u_old=    f(J)+ dt*g(J)+ 1/2*c^2*Tridiag(J+1,s2,-2*s2,s2)*f(J);
u_new=    Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_old  - u_ancient;
c1=0;c2=0;%boundary values
%enforcing the boundary conditions
u_new(1)=c1; u_new(end)=c2;
u_old(1)=c1; u_old(end)=c2;
for i=1:n-2,
%Iterating the system using the Explicit formulation
temp=u_new;
u_new=Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_old - u_ancient;
u_new(1)=c1;  %boundary conditions
u_new(end)=c2;%boundary conditions
u_ancient=u_old;
u_old=temp;
end


That obviously fixes the issue with the boundary points, but I am still getting the wrong picture!! -> http://i.stack.imgur.com/lrVH7.png

It's Working!

I had fixed the issues that was pointed out and the code works fine now! it reads:

    function [u] = Explicit_W(dx,dt,n,initialize_pos,initialize_vel)
%Inputs
%h: length of spatial intervals
X=1;   %X: length of the space domain
J=X/dx; %J: Total # of length intervals
f=str2func(initialize_pos);
g=str2func(initialize_vel);
s=(dt/dx);
c=2;%velocity
r=(s*c)^2;
c1=0;c2=0;%boundary values at the two ends
u=Zero(J+2);u(1)=c1;u(end)=c2;%output solution with its ends fixed
b=Zero(J);b(1)=c1;b(end)=c2;  %vector b enforces the boundary conditions
%u_new and u_old are solution vectors with their ends cut off.
u_old=f(J);%Initializing u(0)
u_new=f(J)+ dt*g(J)+ 1/2*r*(Tridiag(J-1,1,-2,1)*f(J)+b);
for i=1:n-2,
%Iterating the system using the Explicit formulation
temp=u_new;
u_new=r*(Tridiag(J-1,1,-2,1)*u_new+b)+2*eye(J-1)*u_new - u_old;
u_old=temp;
end
u(2:end-1)=u_new;
end

• I think a mod has to take a look at it. I've flagged the question. It shouldn't be long, they're on here often. :) – Michael Brown May 2 '13 at 11:10
• Anyway, it looks like boundary conditions aren't being implemented. What are your boundary conditions? Depending on whether they are Dirichlet or Neumann or something more complicated you have to modify your matrix. There are plenty of references that discuss this, like here for more about the general technique you are using, which is sometimes called the method of lines, and numerical recipes. – Michael Brown May 2 '13 at 11:28
• I thought I have made that clear. I feed in initial values: Dirichlet(initialize_pos, I pass it to f) and a Neumann(initialize_vel, I pass it g). Assigning these two, takes place at the three lines after: %Initializing the system. so I don't know what the problem might be. – Shb May 2 '13 at 11:40
• Those are initial conditions. You also need to impose boundary conditions at every time step. For instance $u(t,0)=1$, $du(t,1)/dx=0$ would be Dirichlet at $x=0$ and Neumann at $x=1$. Your analytic solution is for an unbounded domain (though your numerical solution is on a bounded domain - naturally). If you solved the PDE analytically on the finite domain you would see the need for boundary conditions. – Michael Brown May 2 '13 at 12:11
• oooh! I see... I thought since my initial condition satisfies the boundary conditions given in my problem (namely that u=0 at both ends), that would automatically make the system evolve in a way that would satisfy the boundary conditions. so i guess I should set the first and last elements of u_new explicitly equal to zero after updating it, right? – Shb May 2 '13 at 12:17

There may be several issues at work here, but one problem that stands out is your implementation of boundary conditions.

You're working on the interval $[0,1]$, and partitioning it into $J = 1/dx$ subintervals. Including nodes on the boundary, this gives you $J+1$ total points, which we can call $x_1,\ldots, x_{J+1}$, with $x_1 = 0$ and $x_{J+1} = 1$.

Because you are imposing Dirichlet boundary conditions, the value of $U$ at $x_1$ and $x_{J+1}$ is fixed throughout the simulation, and thus these points should not be included in your update matrix. Your update is of the form

$$U(i+1,:) = (2I + c^2s_2A)U(i,:) - U(i-1,:)$$

where $s_2 = (dt/dx)^2$ and

$$A = \left(\begin{array}{ccccc} -2 & 1 & & & \\ 1 & -2 & 1 & \\ & 1 & -2 & 1\\ &&&\ddots&&\\&&&1&-2\end{array}\right)$$

is the standard Laplacian matrix.

But currently, you are applying this update to the vector $U(i,1:J+1)$, when you should be applying it to $U(i,2:J)$. Put differently, you should only be storing in u_ancient, u_old, and u_new the value of $U$ at $x_2,\ldots,x_{J}$. In fact, this formulation of the problem -- and homogenous Dirichlet boundary conditions -- are encoded directly into the matrix $A$ that you are (implicitly) using. Observe that

\begin{align*} AU(i,2:J) &= \left(\begin{array}{c} -2U(i,1) + U(i,2) \\ U(i,1) - 2U(i,2) + U(i,3) \\ \vdots\\U(i,J-1) - 2U(i,J)\end{array}\right)\\ &=\left(\begin{array}{c} U(i,0) -2U(i,1) + U(i,2) \\ U(i,1) - 2U(i,2) + U(i,3) \\ \vdots\\U(i,J-1) - 2U(i,J) + U(i,J+1)\end{array}\right) \end{align*}

Here the first step is just matrix multiplication, but the second step follows precisely because you have assumed homogenous Dirichlet boundary conditions to begin with, i.e. because $U(i,0) = U(i,J+1) = 0$. This hopefully explains why you needn't include the points $x_1$ and $x_{J+1}$ in your system.

EDIT: Upon re-reading the question, it occurs to me that the issue outlined above shouldn't have altered the solution on account of the problem's explicit formulation, so I'm afraid this doesn't explain your aberrant results.

ANOTHER EDIT: I think I found your problem: the way you've named your variables has led you to implement the time-stepping incorrectly.

Note that after the initial assignment of u_ancient, u_old, and u_new, it is the variable u_new that holds the current value of $U$. In other words u_old is $U(i-1,:)$ and u_new is $U(i,:)$. So in the next step, you want to base the new value, $U(i+1,:)$, off of the current value, $U(i,:)$, and the value that is one step behind, $U(i-1,:)$. You never actually need the value that is two timesteps behind, which you have named u_ancient. In contrast, you've been performing the update

$$U(i+1,:) = (2I+c^2s_2A)U(i-1,:) - U(i-2,:)$$

and this is inconsistent (in the techinical sense) with the PDE.

The fix is very simple. You don't need the variable u_ancient at all, for one. And your loop should be replaced by

for i=1:n-2
temp=u_new
u_new=Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_new - u_old
u_old=temp;
end

• Ha! well spotted sir... this has been driving me nuts in the last 24 hours... :) Thank you so much! – Shb May 2 '13 at 19:05