# Solving PDE in 1D with FD and MATLAB

I have to solve the following equation:

$-u_{xx}=1$,

with $x\in(0,1)$ and $u(0)=0,u(1)=0$. I have to solve it with the following numerical scheme:

$\frac{1}{h_k^2}(-\frac{1}{2}u_{k-1}+u_k-\frac{1}{2}u_{k+1})=1$. So I have to use a non-uniform grid. I have done this so far:

n = 4; %Number points
k = 0:n;
x = 0.5 - 0.5*cos(pi*k/n); %Function to generate points

h = diff(x);
h = h(1:n-1);
h = (h.^2)'; %Difference from point 1 to n-1

b = ones(n-1,1);
b = h.*b; %Solve

A=sparse(diag(2*ones(n-1,1))+diag((-1)*ones(n-2,1),1)+diag((-1)*ones(n-             2,1),-1));

u=A\b;


As you can see I was trying to define everything I need to use a for loop and also trying to do that using just matrix multiplications. The problem is that with for loop I am not able to understand what should be the value at $u(1)$ since I only have $u(0)$ and the method requires three points. With the "matrix version" I am not able to understand how to compute the vector containing the differences between the pints.

• Actually, that's not a PDE but a boundary value problem for an ODE. – Dirk Nov 2 '16 at 22:16

First of all, your finite difference is wrong, as Steve pointed out. It has to be $$\frac{1}{h^2} \left( -u_{i-1} + 2 u_{i} - u_{i+1} \right) = 1.$$ This applies only to a structured grid (with constant space between the points)!

When you want to do the calculation with finite differences, you don't need three points as you dont advance in time. Keep in mind that this is a boundary value problem!

What you need to do is to write down the finite difference equation at every point of your field from 0 to 1. \begin{array} -\frac{1}{h^2} u_0 +\frac{2}{h^2} u_{1} - \frac{1}{h^2} u_{2} = 1 \\ ... \\ -\frac{1}{h^2} u_{n-2} +\frac{2}{h^2} u_{n-1} - \frac{1}{h^2} u_{n} = 1 \end{array} $u_0$ und $u_{n}$ are your boundary conditions.

Now you can create a system of equations like $Au=r$ with $u$ being your solution vector and $r$ a vector containing only ones.

I don't have MATLAB, so I cant help you with the code, but here's an example with Python.

    import numpy as np
import numpy.linalg as LA

n = 21 # number of points
x = np.linspace(0,1,n) # positions of the nodes between 0 and 1
h = (x[-1]-x[0])/(n-1) # distance between two points
A = np.zeros((n-2,n-2))

for iter in range(n-2):
if iter == 0:
A[iter,iter] = 2/h**2
A[iter,iter+1] = -1/h**2
elif iter == n-3:
A[iter,iter] = 2/h**2
A[iter,iter-1] = -1/h**2
else:
A[iter,iter] = 2/h**2
A[iter,iter-1] = -1/h**2
A[iter,iter+1] = -1/h**2

rhs = np.ones((n-2,1))
A_inv = LA.inv(A)
u = np.dot(A_inv,rhs)


Before you use a non-uniform grid, you might try it with a uniform-grid. Don't make it to hard in the beginning.

• Thank you for the reply. The problem is that I have to use a non-homogeneous grid so I was trying to do that from the beginning! Thanks. I did mistake with the derivation now it is correct! – wrong_path Nov 2 '16 at 10:02
• Try to do it with the uniform grid first to understand how to solve the problem in general. To use a non-uniform grid you have to do a transformation of your derivative between your grid and a reference grid, in which you are calculating the derivative. – P. G. Nov 2 '16 at 10:32
• I think I am able to do that with uniform grids (I do not know how to check if the solution is correct), but the problem is only with non-uniform grids. The points are outside the range $0,1$. I am using only n = 4; k = 0:n; x = 1 - 0.5*cos(pi*k/n); as in the notes written by the instructor. – wrong_path Nov 2 '16 at 11:43
• @LorenzoFabbri Your points lie outside $(0,1)$? Well your PDE is defined on $(0,1)$ with BCs on $0$ and $1$, so not having your spatial points in this range doesn't makes sense to me. See my edit answer for a link to non-uniform finite difference operators – Steve Nov 2 '16 at 12:21
• @LorenzoFabbri You can easily check if your solution is correct by using the analytical solution to your problem. If your instructor – P. G. Nov 2 '16 at 18:42

A few points:

• Your points x_k aren't in your domain $(0,1)$. Check plot(x_k): try instead x_k = 0.5 - 0.5*cos(k*pi/n).
• Your vector h could also include the distance to the endpoints 0 and 1
• You can get the distance between points with diff(x_k)
• I think your finite difference operator my be out by a factor of 2. You might want to check that, also what the error is on an irregular mesh
• Check the sizes of all your vectors and matrices: h is 1-by-2, m is 3-by-3 and u is 5-by-1

I have previously answered a question on Finite-difference approximation of the 2nd derivative operator matrix for a staggered grid.

You've said in a comment that you don't know how to check if your solution on a uniform grid is correct: You can investigate this by using the exact solution $$u(x) = \frac{x(1-x)}{2}\,,$$ (found by integrating your ODE twice and imposing the boundary conditions to find the constants of integration), and comparing the numerical solution with the exact solution as $n$ increases, something like $$E_n:= \max_{j\in[1,\dotsc,n]}\left| U_j - u(x_j) \right|\,,$$ where $U_j$ is your numerical solution at the point $x_j$.