The finite difference method for the heat equation-error

[mathmari]

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?

 

[mathmari]

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?

 

[I like Serena]

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]

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]

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?

 

[I like Serena]

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??

 

[mathmari]

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)|$$

 

[I like Serena]

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]

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

 

[I like Serena]

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}$$

 

[I like Serena]

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]

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:

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:

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:

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:

for(j=0; j<J-2; j++){
	UDiag[j]=-(g((j+1)*h+h/2, (n+1)*Dt)/(pow(h,2)));
}

 

[I like Serena]

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:

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]

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?

 

[mathmari]

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?

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;
}

 

[I like Serena]

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]

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!