# Implicit finite difference method

#### mathmari

##### Well-known member
MHB Site Helper
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$?

#### Klaas van Aarsen

##### MHB Seeker
Staff member
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

##### Well-known member
MHB Site Helper
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!!!

#### mathmari

##### Well-known member
MHB Site Helper
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:
Code:
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?

Last edited:

#### mathmari

##### Well-known member
MHB Site Helper
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

#### Klaas van Aarsen

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

#### mathmari

##### Well-known member
MHB Site Helper
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?

#### Klaas van Aarsen

##### MHB Seeker
Staff member
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

##### Well-known member
MHB Site Helper
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}$

#### Klaas van Aarsen

##### MHB Seeker
Staff member
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

##### Well-known member
MHB Site Helper
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$?

#### Klaas van Aarsen

##### MHB Seeker
Staff member
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 $$\displaystyle \frac{1}{Dt^2}+\frac{2}{h^2}$$ and for their neighbors $$\displaystyle - \frac{1}{h^2}$$.
Hmm, you must have rearranged it differently from how I did it.

#### mathmari

##### Well-known member
MHB Site Helper
I rearranged the formula you gave earlier and found for the diagonal elements $$\displaystyle \frac{1}{Dt^2}+\frac{2}{h^2}$$ and for their neighbors $$\displaystyle - \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??

#### Klaas van Aarsen

##### MHB Seeker
Staff member
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??
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

##### Well-known member
MHB Site Helper
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!!!!!

But now having changed some values, the function that calculates $b$ is the following:
Code:
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] ???

Last edited:

#### Klaas van Aarsen

##### MHB Seeker
Staff member
I made a mistake at the way I calculated the vector $b$.
Now I have the same results as you have!!!!!
Good!!!

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

##### Well-known member
MHB Site Helper
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?

Staff member

#### mathmari

##### Well-known member
MHB Site Helper
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.

Ok! Then could I write the function like that?
Code:
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++;
}
}

#### Klaas van Aarsen

##### MHB Seeker
Staff member
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

##### Well-known member
MHB Site Helper
Yep. You can.
But only if the array A has J+1 entries, since you refer to A[J].

At which point do I refer to A[J]??
Do you mean at the following for-loop?
Code:
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]$?

#### Klaas van Aarsen

##### MHB Seeker
Staff member

At which point do I refer to A[J]??
Do you mean at the following for-loop?
Code:
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

##### Well-known member
MHB Site Helper
You're right. So your array A should have J elements.
Does it?
I have defined the array A so:
Code:
double A[J-1];

#### Klaas van Aarsen

##### MHB Seeker
Staff member
I have defined the array A so:
Code:
double A[J-1];
But then A[J-1] is out-of-range!
And it will be unpredictable which results you'll get.

#### mathmari

##### Well-known member
MHB Site Helper
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:
Code:
double A[J];