How to use MATLAB to solve this PDE

Realize first that the equation you have to solve is a linear wave equation, and the numerical diagram you give can be rewritten as

( u^(n+1)_m - 2u^n_m + u^(n-1)_m )/k^2 = ( u^n_(m-1) - 2u^n_m + u^n_(m+1) )/h^2

      

        
        
        
      

    

where k

is the time step, and h

is delta x

in space.

The reformulated numerical scheme clearly shows that the left and right sides are centered second-order finite-difference approximations u_tt

and, u_xx

respectively.

To solve the problem numerically, you need to use the form provided to you, because this is an explicit update formula that you need to implement numerically: it gives you the solution at a point in time n+1

as a function of the previous two times n

and n-1

. You need to start with the initial condition and execute the decision in a timely manner.

Note that the decision imposed on the boundaries of the field ( x=0

and x=1

), so the value of the sampled solutions u^(n)_0

and u^(n)_10

are known for all n

( t=n*k

). At the n

th time step, your unknown vector [u^(n+1)_1, u^(n+1)_2, ..., u^(n+1)_9]

.

Note also that using the update formula to find a solution in a step n+1

requires knowledge of the solution in the previous two steps. So how do you start with n=0

if you need information from the previous two times? This is where initial conditions come into play. You have a solution in n=0

( t=0

), but you also have u_t

at t=0

. These two combinations of information can give you both u^0

and u^1

and get started.

I would use the following startup scheme:

u^0_m = u(h*m,0)  // initial condition on u
(u^2_m - u^0_m)/(2k) = u_t(h*m,0)  // initial condition on u_t

      

        
        
        
      

    

which, combined with the numerical schema used with n=1

, gives you everything you need to define a linear system for u^1_m

and u^2_m

for m=1,...,9

.

Summarizing:

- use the initial scheme to find a solution in n=1

and n=2

at the same time.

- from there on a march in time using the digital circuit you give.

If you are completely lost, check out things like: finite difference schemes, finite difference schemes for advection equations, finite difference schemes for hyperbolic equations, time march.

EDITING:

For a boundary condition on, u_x

you usually use the ghost cell method:

• Imagine a ghost cell at m=-1
  
  , that is, a dummy (or auxiliary) grid point that is used to solve the boundary condition, but this is not part of the solution.

• The first node is m=0
  
  returned to your unknown vector i.e. you are now working with [u_0 u_1 ... u_9]
  
  .

• Use the left boundary condition to close the system. In particular, writing the centered approximation of the boundary condition

  u^n_(1) - u^n_(-1) = 2*h*u_x(0,k*n)
  
  

• The above equation allows you to express a solution to a ghost node in terms of a solution to an internal, real node. So you can apply the digital timing schema (the one you provided) to m=0
  
  node. (The numerical schema applied to m=0
  
  will contain contributions from the m=-1
  
  ghost node, but now you have an expression expressed in terms of m=1
  
  node.)