# The finite difference method for the heat equation-error

#### mathmari

##### Well-known member
MHB Site Helper
Hey!!!

I am implementing in a program the finite difference method for the heat equation.
The problem is the following:
$$u_t(x,t)=(g(x,t)u_x(x,t))_x+f(x,t), \forall (x,t) \in [0,1]x[0,1]$$
$$u(0,t)=u(1,t)=0, \forall t \in [0,1]$$
$$u(x,0)=0, \forall x \in [0,1]$$

where $f(x,t)=\pi x \cos(\pi x t) (1-x)-g(x,t)(2 \pi t \cos(\pi x t)-(x-1) \pi^2 t^2 \sin(\pi x t))+g_x(x,t) (\sin(\pi x t)+\pi (x-1) t \cos(\pi x t))$

$u(x,t)=\sin( \pi x t)(1-x)$

$g(x,t)=2+ \sin(xt)$

The method is the following: ($U_j^n$ is the approximation of $u(x_j,t_n), h=\frac{1}{J}, Dt=\frac{1}{N}, x_j=jh, t_n=nDt$)
$$U_j^0=0, j=0,...,J$$

$$n=0,...,N-1:\frac{U_j^{n+1}-U_j^n}{Dt}=\frac{1}{h}[g(x_j+\frac{h}{2},t_{n+1})\frac{U_{j+1}^{n+1}-U_j^{n+1}}{h}-g(x_j-\frac{h}{2},t_{n+1})\frac{U_j^{n+1}-U_{j-1}^{n+1}}{h}]+f(x_j,t_{n+1}), j=1,...,J-1$$
$$U_0^{n+1}=U_J^{n+1}=0$$

I found the following errors:
For $J=N=10: 0.179505$
For $J=N=20: 0.089506$

Could you tell me if these errors are correct?

Last edited:

#### mathmari

##### Well-known member
MHB Site Helper
I made some modifications at my program, and now I find these errors:
for $J=N=10: 0.174628$
for $J=N=20: 0.086412$

But when we double $J$ and $N$, shouldn't the error get halved?

#### Klaas van Aarsen

##### MHB Seeker
Staff member
Hey!!!

where $f(x,t)=\pi x \cos(\pi x t) (1-x)-g(x,t)(2 \pi t \cos(\pi x t)-(x-1) \pi^2 t^2 \sin(\pi x t))+g_x(x,t) (\sin(\pi x t)+\pi (x-1) t \cos(\pi x t))$
Can it be that it should be:
\begin{aligned}
f(x,t)=\pi x \cos(\pi x t) (1-x)&+\ g(x,t)(2 \pi t \cos(\pi x t)-(x-1) \pi^2 t^2 \sin(\pi x t)) \\
&+\ g_x(x,t) (\sin(\pi x t)+\pi (x-1) t \cos(\pi x t))
\end{aligned}
With a $+$ instead of a $-$?

#### mathmari

##### Well-known member
MHB Site Helper
Hey!!!

Can it be that it should be:
\begin{aligned}
f(x,t)=\pi x \cos(\pi x t) (1-x)&+\ g(x,t)(2 \pi t \cos(\pi x t)-(x-1) \pi^2 t^2 \sin(\pi x t)) \\
&+\ g_x(x,t) (\sin(\pi x t)+\pi (x-1) t \cos(\pi x t))
\end{aligned}
With a $+$ instead of a $-$?
Yes you are right!!!
I replaced the functions at the differential equation and it should be a $+$ instead of a $-$...

#### mathmari

##### Well-known member
MHB Site Helper
After changing the sign of f, I ran again the program and get the following:
for $J=N=4: 0.261844$
for $J=N=8: 0.137994$

Are these errors correct?

#### Klaas van Aarsen

##### MHB Seeker
Staff member
After changing the sign of f, I ran again the program and get the following:
for $J=N=4: 0.261844$
for $J=N=8: 0.137994$

Are these errors correct?
Ha! I think I finally got my expressions right.

With $J=N=4$, I get a maximum deviation in $U_j^N$ of $0.035486$.
With $J=N=8$, that error is $0.00765$.
With $J=N=10$, it is $0.004647$.
With $J=N=20$, it is $0.000768$.

Or were you talking about an other error??

#### mathmari

##### Well-known member
MHB Site Helper
Or were you talking about an other error??
I'm talking about the error that is calculated by $$\underset{1 \leq n \leq N}{max} \underset{1 \leq j \leq J-1}{max}|U_j^n-u(x_j,t_n)|$$

#### Klaas van Aarsen

##### MHB Seeker
Staff member
I'm talking about the error that is calculated by $$\underset{1 \leq n \leq N}{max} \underset{1 \leq j \leq J-1}{max}|U_j^n-u(x_j,t_n)|$$
Oh! Of course.

Anyway, with $J=N=4$, I still get $0.035486$.

#### mathmari

##### Well-known member
MHB Site Helper
Oh! Of course.

Anyway, with $J=N=4$, I still get $0.035486$.
Because my errors differ from yours, I post the approximations that I found for $J=N=5$.. Could you tell if they are also different from yours?

n=0:
0.092009
0.131979
0.119732
0.056693

n=1:
0.185583
0.257637
0.224462
0.106287

n=2:
0.285664
0.383913
0.318851
0.150719

n=3:
0.394673
0.508693
0.395790
0.186965

n=4:
0.509846
0.621373
0.441645
0.208845

Last edited:

#### Klaas van Aarsen

##### MHB Seeker
Staff member
Because my errors differ from yours, I post the approximations that I found for $J=N=5$.. Could you tell if they are also different from yours?
I've got:
$$u = \begin{bmatrix} 0 & 0.100267 & 0.198952 & 0.2945 & 0.385403 & 0.470228 \\ 0 & 0.149214 & 0.289052 & 0.410728 & 0.506597 & 0.570634 \\ 0 & 0.14725 & 0.273819 & 0.361931 & 0.399211 & 0.380423 \\ 0 & 0.096351 & 0.168866 & 0.199605 & 0.180965 & 0.117557 \\ \end{bmatrix}$$
$$U = \begin{bmatrix} 0 & 0.099991 & 0.199442 & 0.298364 & 0.395847 & 0.490025 \\ 0 & 0.148563 & 0.289159 & 0.415125 & 0.518709 & 0.591527 \\ 0 & 0.146334 & 0.273258 & 0.364719 & 0.406867 & 0.390194 \\ 0 & 0.095615 & 0.168179 & 0.2005 & 0.183009 & 0.116612 \\ \end{bmatrix}$$

#### Klaas van Aarsen

##### MHB Seeker
Staff member
Here's what I did:
$$U^{n+1} = (A^{n+1})^{-1}(U^n + \Delta t \cdot f^{n+1})$$

What did you do?

#### mathmari

##### Well-known member
MHB Site Helper
Here's what I did:
$$U^{n+1} = (A^{n+1})^{-1}(U^n + \Delta t \cdot f^{n+1})$$

What did you do?
So you solved a system of the form $Ax=b$, right? I did the same by solving this equation with the Gaussian elimination for tridiagonal systems.
I calculated $b$ as followed:
Code:
for(j=1; j<J-1; j++){
b[j-1]=f(j*h,(n+1)*Dt)+U[SUB]j[/SUB][SUP]n[/SUP]/Dt;
}
But why do we have then different results? Have I maybe done something wrong with the matrix A?

The diagonal of the matrix is:
Code:
for(j=0; j<J-1; j++){
Diag[j]=(g((j+1)*h+h/2, (n+1)*Dt)+g((j+1)*h-h/2, (n+1)*Dt))/(pow(h,2))+1.0/Dt;
}
The elements that are under the main diagonal are:
Code:
LDiag[0]=0;
for(j=1; j<J; j++){
LDiag[j]=-g((j+1)*h-h/2, (n+1)*Dt)/(pow(h,2));
}
And the elements that are above the main diagonal are:
Code:
for(j=0; j<J-2; j++){
UDiag[j]=-(g((j+1)*h+h/2, (n+1)*Dt)/(pow(h,2)));
}

#### Klaas van Aarsen

##### MHB Seeker
Staff member
So you solved a system of the form $Ax=b$, right? I did the same by solving this equation with the Gaussian elimination for tridiagonal systems.
I calculated $b$ as followed:
Code:
for(j=1; j<J-1; j++){
b[j-1]=f(j*h,(n+1)*Dt)+U[SUB]j[/SUB][SUP]n[/SUP]/Dt;
}
But why do we have then different results? Have I maybe done something wrong with the matrix A?
Hold on just there.
You are effectively iterating from j=1 to J-2.
Shouldn't that be from j=1 to J-1?

#### mathmari

##### Well-known member
MHB Site Helper
Hold on just there.
You are effectively iterating from j=1 to J-2.
Shouldn't that be from j=1 to J-1?
With for(j=1; j<J-1; j++) I get for n=0:
0.092009
0.131979
0.119732
0.056693

With for(j=1; j<J; j++) I get for n=0:
0.099991
0.148563
0.146334
0.095615

I calculated it by hand and for n=0 I get(I hope I have done right the calculations):
0.0914
0.1315
0.1270
0.0818

(For $n=0$, $b=f(x_j,t_1)$ since $U^0=0$)

That means that I have made something wrong at the calculation of $b$, or not?

Last edited:

#### mathmari

##### Well-known member
MHB Site Helper
I think that the way that I calculate the function $f$ is wrong, because I get different results as I got when I calculated it by hand..

With the following code I get
3.048142
3.361183
3.327608
2.973833

But it should be
3.00855798
3.30227614071
3.269475746
2.93579504
or not?

Code:
int main(){
int j,J,N;
printf("Give N and J:\n");
scanf("%d %d", &N, &J);

for(j=1; j<J; j++) {
printf("\n %lf\n", f(xx(j,J),tn(1,N)));
}
return 0;
}

double xx(int j,int J){
double x,h=1.0/(double)J;
x=j*h;
return x;
}

double tn(int n, int N){
double t, Dt=1.0/N;
t=n*Dt;
return t;
}

double g(double x, double t){
double gx;
gx=2+sin(x*t);
return gx;
}

double dg(double x, double t){
double dgx;
dgx=t*cos(x*t);
return dgx;
}

double f(double x,double t){
double fx;
double pi=4.0*atan(1.0);
fx=pi*x*cos(pi*x*t)*(1-x)+g(x,t)*(2*pi*t*cos(pi*x*t)-(x-1)*pi*pi*t*t*sin(pi*x*t))+dg(x,t)*(sin(pi*x*t)+pi*(x-1)*t*cos(pi*x*t));

return fx;
}

Last edited:

#### Klaas van Aarsen

##### MHB Seeker
Staff member
With for(j=1; j<J-1; j++) I get for n=0:
0.092009
0.131979
0.119732
0.056693

With for(j=1; j<J; j++) I get for n=0:
0.099991
0.148563
0.146334
0.095615

I calculated it by hand and for n=0 I get(I hope I have done right the calculations):
0.0914
0.1315
0.1270
0.0818

(For $n=0$, $b=f(x_j,t_1)$ since $U^0=0$)

That means that I have made something wrong at the calculation of $b$, or not?
With your correction of $J$ instead of $J-1$, you get the same result as me, so I think that is correct.
With your calculation by hand, I suspect you made a mistake somewhere.

I think that the way that I calculate the function $f$ is wrong, because I get different results as I got when I calculated it by hand..

With the following code I get
3.048142
3.361183
3.327608
2.973833

But it should be
3.00855798
3.30227614071
3.269475746
2.93579504
or not?
Your code looks fine and its result matches mine.
I get:
$$f(x_j, t_1) = \begin{bmatrix} 3.048142 \\ 3.361183 \\ 3.327608 \\ 2.973833 \end{bmatrix}$$

What is your reason to think it should be otherwise?

#### mathmari

##### Well-known member
MHB Site Helper
With your correction of $J$ instead of $J-1$, you get the same result as me, so I think that is correct.
With your calculation by hand, I suspect you made a mistake somewhere.

Your code looks fine and its result matches mine.
I get:
$$f(x_j, t_1) = \begin{bmatrix} 3.048142 \\ 3.361183 \\ 3.327608 \\ 2.973833 \end{bmatrix}$$

What is your reason to think it should be otherwise?
I had only the same results as yours only for $U_j^1$ .
Now I changed some values of the integers at the for loop and I get the same results as yours for all n.

Thank you very much for your help!!!