How to Solve the Time Dependent 1D Heat Equation Using Crank-Nicolson Method?

[apgt512]

Homework Statement

Solve the time dependent 1D heat equation using the Crank-Nicolson method.
The conditions are a interval of length L=1, initial distribution of temperature is u(x,0) = 2-1.5x+sin(pi*x) and the temperature in the ends of the interval are u(0,t) = 2; u(1,t) = 0.5.

Homework Equations

∂u/∂t = √2 (∂²/∂x²) u.

The Attempt at a Solution

I'm really lost here. Since we started seeing finite difference methods and the heat equation I got lost. What I know is: this can be seen as a matrix problem A*u = b. A is a tridiagonal matrix and I have to find the solution matrix u. I think I know what are the elements of the tridiagonal matrix but not sure.
What I don't know is once I get u, that's the temperature at a certain point in time, then how do I solve for the next point in time?. Do I take the u I got as the b for the next step?.
This is the code I've written so far. It uses the LU decomposition method to find the solution matrix. It's written in C.

#include <stdio.h>
#include <math.h>

double * lu(double *diag,double *diag_sup,double *diag_inf,double *b){
    int i;
    double w[n],u[n],y[n];
    static  double x[n];
    w[0]=diag[0];
    for(i=0;i<(n-1);i++){
        u[i]=diag_sup[i]/w[i];
        w[i+1]=diag[i+1]-diag_inf[i+1]*u[i];
    }

    y[0]=b[0]/w[0];    for(i=1;i<n;i++){
        y[i]=(b[i]-y[i-1]*diag_inf[i])/w[i];
    }

    x[n-1]=y[n-1];

    for(i=(n-1);i>=1;i--){
        x[i-1]=y[i-1]-x[i]*u[i-1];
    }

    return x;
}

int main(){
    //Variables
    double alfa,beta,delta,gamma;
    double a[n],e[n],c[n],b[n];
    int i;

    //Boundary conditions
    alfa=0.;
    beta=0.;
    delta=0.;
    gamma=0.;

    //Tridiagonal matrix
    for (i=0; i<n; i++){
        a[i] = 5.;
        e[i] = 2.;
        c[i+1] = 2.;
    }

    //Right side matrix
    for(i=0;i<n;i++){
        b[i]=i;
    }

    //Solution
    double *sol;
    sol = lu(a,e,c,b);
    for(i=0;i<n;i++){
        printf("%lf %lf\n",i*h,*(sol+i));
    }
}

 

[Chestermiller]

Let's see your finite difference equations.

 

[apgt512]

Let's see your finite difference equations.

Well this is the equation:
$$-ru_{i+1}^{n+1} + (1+2r)u_i^{n+1} - ru^{n+1}_{i-1} = ru^n_{i+1} + (1-2r)u_i^n + ru^n_{i-1}.$$
Where
$$r=\frac{\sqrt{2} \Delta t}{2(\Delta x)^2}$$
And that equation can be seen as a matrix problem $Au=b$. But I don't know how to convert the equation to a matrix problem. What do I put in the tridiagonal matrix A, and what do I put in b?. And after I find the solution u, do I use that solution u as the matrix b for the next step in time?.
Thank you!

 

[Orodruin]

I suggest you start with the case of a few discretisation points only (say three) and look at the equations for each point and try to write that equation in matrix form. You also need to find $b$.

And after I find the solution u, do I use that solution u as the matrix b for the next step in time?.

Yes.

 

[Chestermiller]

I checked your difference equation, and it is correct. b is a column vector, and, since each u on the right hand side is already known, b is the entire rhs of the difference equation for the A row centered at i. The dependent variable values at n+1 are the unknowns you are solving for at each time step. (1+2r) runs all the way down the diagonal of the A matrix, and each element along the main diagonal has this value. For each element immediately to the left and right of the main diagonal in a given row of A, the element values are both -r. All the other elements of the A matrix are zero. So you have 3 diagonals running from the upper left corner of A to the lower right corner.

 

[Orodruin]

Also, as a matter of curiosity, note that $b = B u^{n}$ and $Au^{n+1} = b$, where $B$ is another tridiagonal matrix results in $u^{n+1} = A^{-1} B u^n$ and therefore $u^n = (A^{-1}B)^n u^0$.

 

[Chestermiller]

I have a very interesting approach for you to consider that should greatly simplify the numerics, or at least the amount of computation that needs to be done at each time step. I don't know whether this is a standard way of working such a problem.

Divide the finite difference equation by r to obtain:
$$-u_{i+1}^{n+1} + (2+1/r)u_i^{n+1} - u^{n+1}_{i-1} = u^n_{i+1} + (-2+1/r)u_i^n + u^n_{i-1}\tag{1}$$

So now the A matrix has row elements along the three adjacent diagonals of -1, (2+1/r), -1. This is the new A matrix. Eqn. 1 can now be rewritten as:
$$Au^{n+1}=-Au^n+\frac{2}{r}u^n\tag{2}$$
The solution to Eqn. 2 can then be written as: $$u^{n+1}=-u^n+\frac{2}{r}A^{-1}u^n\tag{3}$$
So, according to Eqn. 3, you can get the solution for $u^{n+1}$ by
1. Applying the inverse of the matrix A to $u^n$
2. Multiplying by 2/r
3. Subtracting $u^n$