MHB Implicit finite difference method

[mathmari]

Hey!

I have a implicit finite difference method for the wave equation.

At step 0, we set: $W_j^0=v(x_j), j=0,...,J$

At the step 1, we set: $W_j^1=v(x_j)+Dtu(x_j)+\frac{Dt^2}{2}(\frac{v(x_{j-1})-2v(x_j)+v(x_{j+1})}{h^2}+f(x_j,0)), j=0,...,J$

Can that be that at the step 1 $j$ begins from $0$ and ends at $J$?

 

[I like Serena]

Hey!

I have a implicit finite difference method for the wave equation.

At step 0, we set: $W_j^0=v(x_j), j=0,...,J$

At the step 1, we set: $W_j^1=v(x_j)+Dtu(x_j)+\frac{Dt^2}{2}(\frac{v(x_{j-1})-2v(x_j)+v(x_{j+1})}{h^2}+f(x_j,0)), j=0,...,J$

Can that be that at the step 1 $j$ begins from $0$ and ends at $J$?

Possibly, or perhaps $j=1..J-1$, so that all values are defined.
Otherwise you will need values for $v(x_{-1})$ and $v(x_{J+1})$ that are currently not defined.

 

[mathmari]

Possibly, or perhaps $j=1..J-1$, so that all values are defined.
Otherwise you will need values for $v(x_{-1})$ and $v(x_{J+1})$ that are currently not defined.

Ok! Thank you! (Happy)

 

[mathmari]

I have also an other question..

I want to implement a program in C for this method..

The wave equation is:
$$w_{tt}(x,t)=w_{xx}(x,t)+f(x,t), \forall (x,t) \in [0,1]x[0,1]$$
$$w(0,t)=w(1,t)=0, \forall t \in [0,1]$$
$$w_t(x,0)=u(x), \forall x \in [0,1]$$
$$w(x,0)=v(x), \forall x \in [0,1]$$

where $u(0)=u(1)=0$

The method is:
At step 0, we set: $W_j^0=v(x_j), j=0,...,J$

At the step 1, we set: $$W_j^1=v(x_j)+Dtu(x_j)+\frac{Dt^2}{2}(\frac{v(x_{j-1})-2v(x_j)+v(x_{j+1})}{h^2}+f(x_j,0)), j=1,...,J-1$$

For $n=1,...,N-1$ we do the following:
At the step n+1 we find the $(W_j^{n+1})_{j=1}^{J-1}$ as the solution of the linear system from:
$$\frac{W_j^{n+1}-2W_j^n+W_j^{n-1}}{Dt^2}=\frac{W_{j-1}^{n+1}-2W_j^{n+1}+W_{j+1}^{n+1}}{2h^2}+\frac{W_{j-1}^{n-1}-2W_j^{n-1}+W_{j+1}^{n-1}}{2h^2}+f(x_j,t_n), j=1,..,.J-1$$
$W_0^{n+1}=W_J^{n+1}=0$

The approximations $(W_j^{n-1})_{j=1}^{J-1}$ are saved at the vector $A$, and the $(W_j^{n})_{j=1}^{J-1}$ at $B$.

So we could find the approximations by solving a system of the form $Ax=b$, right?
I calculated the $b$ as followed:

for(j=1; j<J-1; j++){
	b[j-1]=A[j-1]/(2*h*h)-A[j]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+A[j+1]/(2*h*h)+2*B[j-1]/(Dt*Dt)+f(xx(j,J),tn(n,N));
}
b[J-2]=A[J-2]/(2*h*h)-A[J-1]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+2*B[J-2]/(Dt*Dt)+f(xx(J-1,J),tn(n,N));

Is this correct?

$v(x)=0, u(x)= \pi x (1-x)$
$f(x,t)= \pi^2 x^2 (x-1)sin(\pi x t)+2 \pi t cos(\pi x t)+(1-x) \pi^2 t^2 sin( \pi x t)$
the exact solution is $w(x,t)=sin( \pi x t)(1-x)$

Also I want to find the errors ($max_{1 \leq n \leq N} ( max_{1 \leq j \leq J-1} |W_j^n-w(x_j,t_n)|)$).These are what I've found:
For $J=N=10: 0.099211$
For $J=N=20: 0.049875$
Are these errors correct for these J and N?

 

[mathmari]

To check if my program is correct, could you tell me if the following results of the approximations for $J=N=5$ are right?

n=1:
0.198483
0.292070
0.283378
0.181988

n=2:
0.269957
0.424158
0.434060
0.326570

n=3:
0.329534
0.569906
0.628201
0.522819

n=4:
0.397004
0.727531
0.811541
0.764543

 

[I like Serena]

Mmh, for $n=1$ I get:
$$W^1 = \begin{bmatrix}
0 \\
0.100531 \\
0.150796 \\
0.150796 \\
0.100531 \\
0
\end{bmatrix}$$

 

[mathmari]

Mmh, for $n=1$ I get:
$$W^1 = \begin{bmatrix}
0 \\
0.100531 \\
0.150796 \\
0.150796 \\
0.100531 \\
0
\end{bmatrix}$$

Ok.. But when we solve the system do we not get as answer the approximations $W^{n+1}$? So for $n=1$ we get $W^2$, or did I misunderstand it?
This $W^1$ is from the step 1, right?

 

[I like Serena]

Ok.. But when we solve the system do we not get as answer the approximations $W^{n+1}$? So for $n=1$ we get $W^2$, or did I misunderstand it?
This $W^1$ is from the step 1, right?

Yep.

For more steps, I get:
$$W_{j=1..J-1}^n = \begin{bmatrix}
0 & 0.100531 & 0.167758 & 0.190311 & 0.183852 & 0.178421 \\
0 & 0.150796 & 0.252344 & 0.279416 & 0.244453 & 0.191595 \\
0 & 0.150796 & 0.246062 & 0.256217 & 0.191528 & 0.095445 \\
0 & 0.100531 & 0.156116 & 0.146507 & 0.08322 & -0.002342 \\
\end{bmatrix}$$

Hmm... that is not the same as what you have...

 

[mathmari]

Yep.

For more steps, I get:
$$W_{j=1..J-1}^n = \begin{bmatrix}
0 & 0.100531 & 0.167758 & 0.190311 & 0.183852 & 0.178421 \\
0 & 0.150796 & 0.252344 & 0.279416 & 0.244453 & 0.191595 \\
0 & 0.150796 & 0.246062 & 0.256217 & 0.191528 & 0.095445 \\
0 & 0.100531 & 0.156116 & 0.146507 & 0.08322 & -0.002342 \\
\end{bmatrix}$$

Hmm... that is not the same as what you have...

I calculated the approximation $W_j^2$ by hand and I get the following system ( I hope I made no mistakes at the calculations ):
$\begin{bmatrix}
50 & -12.5 & 0 & 0\\
-12.5 &50 &-12.5 &0 \\
0 & -12.5 & 50 &-12.5 \\
0 &0 & -12.5 & 50
\end{bmatrix} W_j^2=\begin{bmatrix}
6.271728\\
8.580236 \\
8.243136 \\
5.557183
\end{bmatrix}$

So $W_j^2=\begin{bmatrix}
0.1984\\
0.2921 \\
0.2834 \\
0.1820
\end{bmatrix}$

 

[I like Serena]

I calculated the approximation $W_j^2$ by hand and I get the following system ( I hope I made no mistakes at the calculations ):
$\begin{bmatrix}
50 & -12.5 & 0 & 0\\
-12.5 &50 &-12.5 &0 \\
0 & -12.5 & 50 &-12.5 \\
0 &0 & -12.5 & 50
\end{bmatrix} W_j^2=\begin{bmatrix}
6.271728\\
8.580236 \\
8.243136 \\
5.557183
\end{bmatrix}$

So $W_j^2=\begin{bmatrix}
0.1984\\
0.2921 \\
0.2834 \\
0.1820
\end{bmatrix}$

I've got:
$$\begin{bmatrix}
75 & -25 & 0 & 0\\
-25 &75 &-25 &0 \\
0 & -25 & 75 &-25 \\
0 &0 & -25 & 75
\end{bmatrix} W_j^2=\begin{bmatrix}
6.271728\\
8.580236 \\
8.243136 \\
5.557183
\end{bmatrix}$$

 

[mathmari]

I've got:
$$\begin{bmatrix}
75 & -25 & 0 & 0\\
-25 &75 &-25 &0 \\
0 & -25 & 75 &-25 \\
0 &0 & -25 & 75
\end{bmatrix} W_j^2=\begin{bmatrix}
6.271728\\
8.580236 \\
8.243136 \\
5.557183
\end{bmatrix}$$

The elements of the diagonal of the matrix are the coefficients of $W_j^{n+1}$, aren't they? So they are equal to $\frac{1}{Dt^2}+\frac{1}{h^2}$.
The elements of the upper diagonal are the coefficients of $W_{j+1}^{n+1}$, so they are equal to $- \frac{1}{2h^2}$.
And the elements of the lower diagonal are the coefficients of $W_{j-1}^{n+1}$, so they are equal to $- \frac{1}{2h^2}$.

Have I calculated something wrong?

Could you tell me how you found that the elements of the diagonal are $75$ and the elements of the lower and upper diagonal $-25$?

 

[I like Serena]

The elements of the diagonal of the matrix are the coefficients of $W_j^{n+1}$, aren't they? So they are equal to $\frac{1}{Dt^2}+\frac{1}{h^2}$.
The elements of the upper diagonal are the coefficients of $W_{j+1}^{n+1}$, so they are equal to $- \frac{1}{2h^2}$.
And the elements of the lower diagonal are the coefficients of $W_{j-1}^{n+1}$, so they are equal to $- \frac{1}{2h^2}$.
Have I calculated something wrong?

Could you tell me how you found that the elements of the diagonal are $75$ and the elements of the lower and upper diagonal $-25$?

I rearranged the formula you gave earlier and found for the diagonal elements $$\frac{1}{Dt^2}+\frac{2}{h^2}$$ and for their neighbors $$- \frac{1}{h^2}$$.
Hmm, you must have rearranged it differently from how I did it.

 

[mathmari]

I rearranged the formula you gave earlier and found for the diagonal elements $$\frac{1}{Dt^2}+\frac{2}{h^2}$$ and for their neighbors $$- \frac{1}{h^2}$$.
Hmm, you must have rearranged it differently from how I did it.

From the following formula, do we not use the terms in red?

$\frac{W_j^{n+1}}{Dt^2}$$+\frac{-2W_j^n+W_j^{n-1}}{Dt^2}=$$\frac{W_{j-1}^{n+1}}{2h^2} - \frac{2W_j^{n+1}}{2h^2} + \frac{W_{j+1}^{n+1}}{2h^2}$$+\frac{W_{j-1}^{n-1}-2W_j^{n-1}+W_{j+1}^{n-1}}{2h^2}+f(x_j,t_n)$

So $\frac{W_j^{n+1}}{Dt^2}-\frac{W_{j-1}^{n+1}}{2h^2}+ \frac{W_j^{n+1}}{h^2}- \frac{W_{j+1}^{n+1}}{2h^2}$$=\frac{2W_j^n-W_j^{n-1}}{Dt^2}+\frac{W_{j-1}^{n-1}-2W_j^{n-1}+W_{j+1}^{n-1}}{2h^2}+f(x_j,t_n)$

$W_j^{n+1}(\frac{1}{Dt^2}+\frac{1}{h^2})+W_{j-1}^{n+1}(\frac{-1}{2h^2})+ W_{j+1}^{n+1}(\frac{-1}{2h^2})=\frac{2W_j^n-W_j^{n-1}}{Dt^2}+\frac{W_{j-1}^{n-1}-2W_j^{n-1}+W_{j+1}^{n-1}}{2h^2}+f(x_j,t_n)$

How did you rearrange it??

 

[I like Serena]

From the following formula, do we not use the terms in red?

$\frac{W_j^{n+1}}{Dt^2}$$+\frac{-2W_j^n+W_j^{n-1}}{Dt^2}=$$\frac{W_{j-1}^{n+1}}{2h^2} - \frac{2W_j^{n+1}}{2h^2} + \frac{W_{j+1}^{n+1}}{2h^2}$$+\frac{W_{j-1}^{n-1}-2W_j^{n-1}+W_{j+1}^{n-1}}{2h^2}+f(x_j,t_n)$

So $\frac{W_j^{n+1}}{Dt^2}-\frac{W_{j-1}^{n+1}}{2h^2}+ \frac{W_j^{n+1}}{h^2}- \frac{W_{j+1}^{n+1}}{2h^2}$$=\frac{2W_j^n-W_j^{n-1}}{Dt^2}+\frac{W_{j-1}^{n-1}-2W_j^{n-1}+W_{j+1}^{n-1}}{2h^2}+f(x_j,t_n)$

$W_j^{n+1}(\frac{1}{Dt^2}+\frac{1}{h^2})+W_{j-1}^{n+1}(\frac{-1}{2h^2})+ W_{j+1}^{n+1}(\frac{-1}{2h^2})=\frac{2W_j^n-W_j^{n-1}}{Dt^2}+\frac{W_{j-1}^{n-1}-2W_j^{n-1}+W_{j+1}^{n-1}}{2h^2}+f(x_j,t_n)$

How did you rearrange it??

I made a mistake.
Now I get the same matrix as you have.

As a result I get:
$$W = \begin{bmatrix}
0 & 0.100531 & 0.198483 & 0.293286 & 0.386943 & 0.482468 \\
0 & 0.150796 & 0.29207 & 0.416943 & 0.521666 & 0.603338 \\
0 & 0.150796 & 0.283378 & 0.38187 & 0.433659 & 0.42937 \\
0 & 0.100531 & 0.181988 & 0.225558 & 0.216531 & 0.15098 \\
\end{bmatrix}$$

... but that is still not quite the same as what you had...

 

[mathmari]

I made a mistake.
Now I get the same matrix as you have.

As a result I get:
$$W = \begin{bmatrix}
0 & 0.100531 & 0.198483 & 0.293286 & 0.386943 & 0.482468 \\
0 & 0.150796 & 0.29207 & 0.416943 & 0.521666 & 0.603338 \\
0 & 0.150796 & 0.283378 & 0.38187 & 0.433659 & 0.42937 \\
0 & 0.100531 & 0.181988 & 0.225558 & 0.216531 & 0.15098 \\
\end{bmatrix}$$

... but that is still not quite the same as what you had...

I made a mistake at the way I calculated the vector $b$.
Now I have the same results as you have! (Smirk)

(Wait)
But now having changed some values, the function that calculates $b$ is the following:

void *b(double A[],double B[],double C[],int J, int N, int n){
	double *pointer=&C[0],Dt=T/(double)N, h=(c-a)/(double)J;
	int j;	
	for(j=1; j<J; j++){
		C[j-1]=A[j-2]/(2*h*h)-A[j-1]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+A[j]/(2*h*h)+2*B[j-1]/(Dt*Dt)+f(xx(j,J),tn(n,N));
		pointer=&C[j-1];
		//printf("%lf\n",*pointer);
		pointer++;		
	}
}

How can that be that it works, when at the first iteration for j=1 we have A[j-2]=A[-1] ?

 

[I like Serena]

I made a mistake at the way I calculated the vector $b$.
Now I have the same results as you have! (Smirk)

Good!

(Wait)
How can that be that it works, when at the first iteration for j=1 we have A[j-2]=A[-1] ?

Aren't all values of A zero in the first iteration?
Then it won't matter much which index on A you are using.
However, the value for A[-1] is undefined and will often contain an unpredictable value.
It seems you've been in "luck" as it would appear A[-1] turned out to be zero.

 

[mathmari]

Aren't all values of A zero in the first iteration?
Then it won't matter much which index on A you are using.
However, the value for A[-1] is undefined and will often contain an unpredictable value.
It seems you've been in "luck" as it would appear A[-1] turned out to be zero.

So should I change these indices or let the for loop be so as it is now?

 

[I like Serena]

So should I change these indices or let the for loop be so as it is now?

Slightly modified quote from The Ten Commandments for C Programmers:

Thou shalt not index an array outside of its specified range, for chaos and madness await thee.​

(Evilgrin) (Swearing) (Angel)

 

[mathmari]

Slightly modified quote from The Ten Commandments for C Programmers:

Thou shalt not index an array outside of its specified range, for chaos and madness await thee.​

(Evilgrin) (Swearing) (Angel)

Ok! Then could I write the function like that?

void *b(double A[],double B[],double C[],int J, int N, int n){
	double *pointer=&C[0],Dt=T/(double)N, h=(c-a)/(double)J;
	int j;	
	C[0]=-A[0]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+A[1]/(2*h*h)+2*B[0]/(Dt*Dt)+f(xx(1,J),tn(n,N));
	pointer=&C[0];
	printf("%lf\n",*pointer);
	for(j=2; j<J; j++){
		C[j-1]=A[j-2]/(2*h*h)-A[j-1]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+A[j]/(2*h*h)+2*B[j-1]/(Dt*Dt)+f(xx(j,J),tn(n,N));
		pointer=&C[j-1];
		printf("%lf\n",*pointer);
		pointer++;		
	}
}

 

[I like Serena]

Ok! Then could I write the function like that?

Yep. You can.
But only if the array A has J+1 entries, since you refer to A[J].

 

[mathmari]

Yep. You can.
But only if the array A has J+1 entries, since you refer to A[J].

(Thinking)
At which point do I refer to A[J]??
Do you mean at the following for-loop?

for(j=2; j<J; j++)
	C[j-1]=A[j-2]/(2*h*h)-A[j-1]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+A[j]/(2*h*h)+2*B[j-1]/(Dt*Dt)+f(xx(j,J),tn(n,N));

But isn't it $j=2,...,J-1$, so $A[J-1]$?

 

[I like Serena]

(Thinking)
At which point do I refer to A[J]??
Do you mean at the following for-loop?

for(j=2; j<J; j++)
	C[j-1]=A[j-2]/(2*h*h)-A[j-1]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+A[j]/(2*h*h)+2*B[j-1]/(Dt*Dt)+f(xx(j,J),tn(n,N));

But isn't it $j=2,...,J-1$, so $A[J-1]$?

You're right. So your array A should have J elements.
Does it?

 

[mathmari]

You're right. So your array A should have J elements.
Does it?

I have defined the array A so:

double A[J-1];

 

[I like Serena]

I have defined the array A so:

double A[J-1];

But then A[J-1] is out-of-range!
And it will be unpredictable which results you'll get.

 

[mathmari]

But then A[J-1] is out-of-range!
And it will be unpredictable which results you'll get.

A ok..So you mean that I should define the array so:

double A[J];

 

[I like Serena]

A ok..So you mean that I should define the array so:

double A[J];

Yes.
And you should also make sure it is initialized to zero.
Is it?

 

[mathmari]

Yes.
And you should also make sure it is initialized to zero.
Is it?

I have initialized it as followed:

for(j=0; j<J; j++){
	A[j]=0;
}

 

[I like Serena]

I have initialized it as followed:

for(j=0; j<J; j++){
	A[j]=0;
}

Good! ;)

 

[mathmari]

Good! ;)

I read again the exercise and I saw that we need an array $A$ with $J-1$ elements since the zero elements don't have to be saved. That means that I have to define this array so: double A[J-1]; , or not?

So I changed the function that calculates $b$:

void *b(double A[],double B[],double C[],int J, int N, int n){
	double *pointer=&C[0],Dt=T/(double)N, h=(c-a)/(double)J;
	int j;	
	C[0]=-A[0]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+A[1]/(2*h*h)+2*B[0]/(Dt*Dt)+f(xx(1,J),tn(n,N));
	pointer=&C[0];
	printf("%lf\n",*pointer);
	for(j=2; j<J-1; j++){
		C[j-1]=A[j-2]/(2*h*h)-A[j-1]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+A[j]/(2*h*h)+2*B[j-1]/(Dt*Dt)+f(xx(j,J),tn(n,N));
		pointer=&C[j-1];
		printf("%lf\n",*pointer);
		pointer++;		
	}
	C[J-2]=A[J-3]/(2*h*h)-A[J-2]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+2*B[J-2]/(Dt*Dt)+f(xx(J-1,J),tn(n,N));
	pointer=&C[J-2];
	printf("%lf\n",*pointer);
}

Now there is no $A[J-1]$..Is it correct now?

 

[I like Serena]

I read again the exercise and I saw that we need an array $A$ with $J-1$ elements since the zero elements don't have to be saved. That means that I have to define this array so: double A[J-1]; , or not?

So I changed the function that calculates $b$:

void *b(double A[],double B[],double C[],int J, int N, int n){
	double *pointer=&C[0],Dt=T/(double)N, h=(c-a)/(double)J;
	int j;	
	C[0]=-A[0]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+A[1]/(2*h*h)+2*B[0]/(Dt*Dt)+f(xx(1,J),tn(n,N));
	pointer=&C[0];
	printf("%lf\n",*pointer);
	for(j=2; j<J-1; j++){
		C[j-1]=A[j-2]/(2*h*h)-A[j-1]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+A[j]/(2*h*h)+2*B[j-1]/(Dt*Dt)+f(xx(j,J),tn(n,N));
		pointer=&C[j-1];
		printf("%lf\n",*pointer);
		pointer++;		
	}
	C[J-2]=A[J-3]/(2*h*h)-A[J-2]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+2*B[J-2]/(Dt*Dt)+f(xx(J-1,J),tn(n,N));
	pointer=&C[J-2];
	printf("%lf\n",*pointer);
}

Now there is no $A[J-1]$..Is it correct now?

Yep. That is correct.

Btw, any reason that the return value is "void *" instead of "void"?
And that you execute "pointer++" when you make no use of the incremented value?

 

[mathmari]

Yep. That is correct.

Great! (Yes) Thank you for your help! (Happy)

Btw, any reason that the return value is "void *" instead of "void"?
And that you execute "pointer++" when you make no use of the incremented value?

I used pointers to return the array and use it in the main function..