# Time dependent heat equation

Tags:
1. Nov 20, 2016

### apgt512

1. The problem statement, all variables and given/known data
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.

2. Relevant equations
∂u/∂t = √2 (∂2/∂x2) u.

3. 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.
Code (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));
}
}

2. Nov 21, 2016

### Staff: Mentor

Let's see your finite difference equations.

3. Nov 21, 2016

### apgt512

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!

4. Nov 21, 2016

### Orodruin

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

Yes.

5. Nov 21, 2016

### Staff: Mentor

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.

6. Nov 21, 2016

### Orodruin

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

7. Nov 22, 2016

### Staff: Mentor

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$