How to use MATLAB to solve this PDE

I have the following question in my exam:

enter image description here

I need to use MATLAB to solve it. The problem is I haven't seen such an issue before and I am struggling to get started.

I have a 1x1 grid split into 10x10. I know that I can calculate the entire bottom row except corners using 1/10 * x * 2. I also know that I can calculate the entire right row using (1/10) (1 + t) ^ 2. However, I can't figure out how to get enough points to be able to fill in values ​​for the entire grid. I know this must have something to do with the partial derivatives given in the problem, but I'm not entirely sure where they come into play (especially the u_x equation). Can someone help me get started here?

I don't need the whole solution. Once I have enough points, I can easily write a Matlab program to do the rest. Actually, I think I just want the x = 0 lattice, then I'll just fill in the middle of the grid.

I calculated the bottom row, minus the two angles, to be 0.001, 0.004, 0.009, 0.016, 0.025, 0.036, 0.049, 0.064, 0.081. And similarly, the entire right-hand line is trivial to compute using the given boundary condition. I just can't get together where to go from there.

Edit: the third boundary state equation was wrong. it should read:

u_x (0, t) = 1 / 5t, NOT u (0, t) = 1 / 5t


source to share

1 answer

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


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



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


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




All Articles