Welcome to our community

Be a part of something great, join today!

Implicit finite difference method

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
Hey!! :eek:

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
Mar 5, 2012
8,780
Hey!! :eek:

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
Apr 14, 2013
4,047
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

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
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
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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... :eek:
 

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
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... :eek:
I calculated the approximation $W_j^2$ by hand and I get the following system ( I hope I made no mistakes at the calculations :eek: ):
$\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
Mar 5, 2012
8,780
I calculated the approximation $W_j^2$ by hand and I get the following system ( I hope I made no mistakes at the calculations :eek: ):
$\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}$$

:eek:
 

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
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}$$

:eek:
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
Mar 5, 2012
8,780
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. :eek:
 

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
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. :eek:
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
Mar 5, 2012
8,780
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. :eek:
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... :eek:
 

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
I made a mistake. :eek:
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... :eek:
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:
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] ??? :confused:
 
Last edited:

Klaas van Aarsen

MHB Seeker
Staff member
Mar 5, 2012
8,780
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] ??? :confused:
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
Apr 14, 2013
4,047
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?
 

Klaas van Aarsen

MHB Seeker
Staff member
Mar 5, 2012
8,780

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
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?
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
Mar 5, 2012
8,780
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]. :eek:
 

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
Yep. You can.
But only if the array A has J+1 entries, since you refer to A[J]. :eek:
(Thinking)
At which point do I refer to A[J]?? :confused:
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
Mar 5, 2012
8,780
(Thinking)
At which point do I refer to A[J]?? :confused:
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
Apr 14, 2013
4,047

Klaas van Aarsen

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

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
But then A[J-1] is out-of-range! :eek:
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];