The finite difference method for the heat equation-error

Click For Summary

Discussion Overview

The discussion revolves around the implementation of the finite difference method for solving the heat equation, specifically addressing errors encountered in numerical approximations. Participants explore the formulation of the equation, modifications to the function definitions, and the resulting impact on error calculations. The scope includes technical reasoning and mathematical formulation.

Discussion Character

  • Technical explanation
  • Mathematical reasoning
  • Debate/contested

Main Points Raised

  • Initial errors reported for different grid sizes are questioned regarding their correctness and expected behavior when grid sizes are doubled.
  • Some participants propose a correction to the formulation of the function \( f(x,t) \) by suggesting a change from a minus to a plus sign in its definition.
  • After modifying the sign in \( f \), new error values are presented, leading to further inquiries about their correctness.
  • Participants share their computed values for different grid sizes and express confusion over discrepancies in results, prompting discussions about the calculation of the matrix \( A \) and the vector \( b \) in the numerical method.
  • There is a suggestion that the iteration limits in the code may have been incorrectly set, impacting the results.
  • Concerns are raised about the accuracy of the function \( f \) based on hand calculations versus programmatic outputs, leading to further scrutiny of the implementation details.

Areas of Agreement / Disagreement

Participants express uncertainty about the correctness of their error calculations and the formulation of the function \( f(x,t) \). There is no consensus on the correct values or the impact of the proposed changes, indicating multiple competing views and unresolved issues.

Contextual Notes

Discrepancies in error calculations and function definitions highlight potential misunderstandings in the implementation of the finite difference method. The discussion reflects a reliance on specific coding practices and mathematical formulations that may not be universally agreed upon.

mathmari
Gold Member
MHB
Messages
4,984
Reaction score
7
Hey! :o

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:
Physics news on Phys.org
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?
 
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 $-$?
 
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! :o
I replaced the functions at the differential equation and it should be a $+$ instead of a $-$...
 
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?
 
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:
 
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)|$$
 
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. :o

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

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? :o
 
  • #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? :o

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

Similar threads

Undergrad Numerically solving a non-local PDE
  • · Replies 3 ·
Replies
3
Views
3K
Undergrad On vibrating string and differentiating infinite sum
  • · Replies 7 ·
Replies
7
Views
3K
Undergrad 1D time dependent PDE - using orthogonal collocation
  • · Replies 2 ·
Replies
2
Views
4K
Graduate Numerically solving a transport equation
  • · Replies 3 ·
Replies
3
Views
3K
Undergrad Derivation of Lax-Wendroff finite volume scheme
  • · Replies 1 ·
Replies
1
Views
3K
Undergrad On the approximate solution obtained through Euler's method
  • · Replies 1 ·
Replies
1
Views
3K
Undergrad Problem getting the coefficients of a non-homogeneous PDE using the Fourier method
  • · Replies 1 ·
Replies
1
Views
2K
MHB Approximating Equations for Unknowns: How to Justify the Form of $U_0$?
  • · Replies 5 ·
Replies
5
Views
3K
MHB Explore the Multiple Grid Method: Get Accurate Solution Faster
  • · Replies 36 ·
2
Replies
36
Views
8K
Undergrad Finite difference scheme for u'=u^2, u0=1, 0<t<2
  • · Replies 8 ·
Replies
8
Views
4K