MHB The finite difference method for the heat equation-error

AI Thread Summary
The discussion revolves around implementing the finite difference method for the heat equation and addressing errors in the numerical solution. Initially, the user encountered discrepancies in error values when doubling the grid sizes, leading to a realization that a sign in the function f was incorrect. After correcting this, the user reported improved error values but sought validation of these results. Further discussions revealed issues in the calculation of the vector b used in the system of equations, with participants comparing their outputs and debugging their code. Ultimately, the user achieved consistent results with the help of the community, confirming the accuracy of their implementation.
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:
Mathematics 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
 
Back
Top