Solving PDE in 1D with FD and MATLAB

Question

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.

Answer 1

First of all, your finite difference is wrong, as Steve pointed out. It has to be \begin{equation} \frac{1}{h^2} \left( -u_{i-1} + 2 u_{i} - u_{i+1} \right) = 1. \end{equation} 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.

Answer 2

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