Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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)
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;
%Initializing the system, 
%u_ancient, u_old and u_new are U(i-1,:),U(i-1,:) and U(i+1,:) respectively.
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
    u_new=Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_old - u_ancient;


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

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)
for j=1:(L+1),


function [ u0 ] = Zero( L)
%Summary of this function goes here
%   Initializes the array to a zero!


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


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


I get funny results, ranging from -124 to 124! Here is a picture of the results:

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?


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_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
    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

That obviously fixes the issue with the boundary points, but I am still getting the wrong picture!! ->

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)
%h: length of spatial intervals
X=1;   %X: length of the space domain
J=X/dx; %J: Total # of length intervals
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
    u_new=r*(Tridiag(J-1,1,-2,1)*u_new+b)+2*eye(J-1)*u_new - u_old;
share|improve this question

migrated from May 2 '13 at 12:44

This question came from our site for active researchers, academics and students of physics.

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
up vote 3 down vote accepted

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
    u_new=Tridiag(J+1,c^2*s2,2-2*c^2*s2,c^2*s2)*u_new - u_old
share|improve this answer
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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.