The finite difference method for the heat equation-error

In summary: With for(j=1; j<J-1; j++) I get for n=0:0.0920090.1319790.1197320.056693With for(j=1; j<J; j++) I get for n=0:0.0999910.1485630.1463340.095615I calculated it by hand and for n=0 I get(I hope I have done right the calculations):0.09140.13150.12700.0818(For $n=0$, $b=f(x_j,t_1)$ since $U^0=0$)That means that I have
  • #1
mathmari
Gold Member
MHB
5,049
7
Hey! :eek:

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 by a moderator:
Mathematics news on Phys.org
  • #2
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?
 
  • #3
Hey! :)

mathmari said:
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 $-$?
 
  • #4
I like Serena said:
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! :eek:
I replaced the functions at the differential equation and it should be a $+$ instead of a $-$...
 
  • #5
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?
 
  • #6
mathmari said:
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. (Cool)

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?? :confused:
 
  • #7
I like Serena said:
Or were you talking about an other error?? :confused:

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)|$$
 
  • #8
mathmari said:
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. :eek:

Anyway, with $J=N=4$, I still get $0.035486$.
 
  • #9
I like Serena said:
Oh! Of course. :eek:

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 by a moderator:
  • #10
mathmari said:
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}$$
 
  • #11
Here's what I did:
$$U^{n+1} = (A^{n+1})^{-1}(U^n + \Delta t \cdot f^{n+1})$$

What did you do? :eek:
 
  • #12
I like Serena said:
Here's what I did:
$$U^{n+1} = (A^{n+1})^{-1}(U^n + \Delta t \cdot f^{n+1})$$

What did you do? :eek:

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)));
}
 
  • #13
mathmari said:
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?
 
  • #14
I like Serena said:
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 by a moderator:
  • #15
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 by a moderator:
  • #16
mathmari said:
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. :eek:
mathmari said:
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?
 
  • #17
I like Serena said:
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. :eek:

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! :eek:
 

FAQ: The finite difference method for the heat equation-error

What is the finite difference method?

The finite difference method is a mathematical technique used to approximate solutions to differential equations. It involves dividing a continuous domain into a discrete grid and using finite differences to approximate the derivatives in the equation.

How does the finite difference method solve the heat equation?

The finite difference method solves the heat equation by discretizing the domain into a grid of points and approximating the derivatives in the equation using finite differences. This results in a system of linear equations that can be solved numerically to obtain an approximate solution to the heat equation.

What is the error in the finite difference method for the heat equation?

The error in the finite difference method for the heat equation refers to the difference between the exact solution of the heat equation and the approximate solution obtained using the finite difference method. This error can be reduced by using a finer grid or by using higher-order finite difference schemes.

What are the advantages of using the finite difference method for the heat equation?

One advantage of using the finite difference method for the heat equation is that it is relatively easy to implement and can handle complex boundary conditions. It is also a versatile method that can be used for a wide range of differential equations. Additionally, the finite difference method is computationally efficient and can handle non-uniform grids.

What are the limitations of the finite difference method for the heat equation?

One limitation of the finite difference method for the heat equation is that it can only approximate solutions on a discrete grid, which may not accurately capture the behavior of the system in between grid points. Additionally, the accuracy of the finite difference method depends on the choice of grid spacing, and using a finer grid can increase the computational cost. It is also limited to problems with simple geometries and may struggle with highly non-linear equations.

Similar threads

MHB Approximating Equations for Unknowns: How to Justify the Form of $U_0$?
Replies
5
Views
3K
MHB Explore the Multiple Grid Method: Get Accurate Solution Faster
2
Replies
36
Views
7K
MHB How is the Finite Difference Method Applied in Wave Approximation?
Replies
1
Views
2K
MHB How to Solve a Boundary Value Problem Using the Finite Difference Method?
Replies
1
Views
2K
MHB Order of Accuracy for Finite Difference Method Backward Euler
Replies
1
Views
2K
MHB Implicit finite difference method
Replies
30
Views
6K
MHB Bounds of the difference of a bounded band-limited function
Replies
1
Views
2K
MHB Interval with Dirac function in a finite interval
Replies
2
Views
1K
MHB How to get the desired upper bound
Replies
1
Views
2K
I Alternating Harmonic Numbers are cool, spread the word!
Replies
4
Views
1K

Hot Threads

Recent Insights

Back
Top