# Finite difference method-convergence

#### mathmari

##### Well-known member
MHB Site Helper
Hello!
I am implementing the finite difference method in a program in C and I got stuck at the rate of convergence.. The formula is $$\frac{log(\frac{e_{1}}{e_{2}})}{log(\frac{J_{2}}{J_{1}})}$$, right? where $$e_{i}=max|y_{j}^{J_{i}}-y(x_{j}^{J_{i}})|$$, $$0<=j<=J_{i}$$. How can I find the $$J_{1}$$ and $$J_{2}$$?

Last edited:

#### mathmari

##### Well-known member
MHB Site Helper
I would also know how the error changes while the step size is getting smaller? Is it getting smaller or bigger??

#### Klaas van Aarsen

##### MHB Seeker
Staff member
Hello!
I am implementing the finite difference method in a program in C and I got stuck at the rate of convergence.. The formula is $$\frac{log(\frac{e_{1}}{e_{2}})}{log(\frac{J_{2}}{J_{1}})}$$, right? where $$e_{i}=max|y_{j}^{J_{i}}-y(x_{j}^{J_{i}})|$$, $$0<=j<=J_{i}$$. How can I find the $$J_{1}$$ and $$J_{2}$$?
Which finite difference method do you mean?
I looked it up on wiki and found a family of methods.
In particular wiki does not refer to $J_i$ or $y_{j}^{J_{i}}$.

#### mathmari

##### Well-known member
MHB Site Helper
Which finite difference method do you mean?
I looked it up on wiki and found a family of methods.
In particular wiki does not refer to $J_i$ or $y_{j}^{J_{i}}$.
I want to implement a finite difference method for a forth order problem.
$$y'''-q(x)y''(x)=f(x), \forall x \in [0,1]$$
$$y(0)=y(1)=y''(0)=y''(1)=0$$
where $f:[0,1] \rightarrow R$, $q:[0,1] \rightarrow R$, $q(x) \geq 0, \forall x \in [0,1]$

To find the solution y, we consider the function $w:[0,1] \rightarrow R$, where $w(x)=y''(x)$, so we have the following problem:
$w''(x)-q(x)w(x)=f(x), \forall x \in [0,1]$ (1)
$w(0)=w(1)=0$
and
$y''=w(x), \forall x \in [0,1]$
$y(0)=y(1)=0$

For the problem (1) when we halve the step, the error get divided by 4, right?

But what happens for the whole problem??

#### Klaas van Aarsen

##### MHB Seeker
Staff member
I want to implement a finite difference method for a forth order problem.
$$y'''-q(x)y''(x)=f(x), \forall x \in [0,1]$$
$$y(0)=y(1)=y''(0)=y''(1)=0$$
where $f:[0,1] \rightarrow R$, $q:[0,1] \rightarrow R$, $q(x) \geq 0, \forall x \in [0,1]$

To find the solution y, we consider the function $w:[0,1] \rightarrow R$, where $w(x)=y''(x)$, so we have the following problem:
$w''(x)-q(x)w(x)=f(x), \forall x \in [0,1]$ (1)
$w(0)=w(1)=0$
and
$y''=w(x), \forall x \in [0,1]$
$y(0)=y(1)=0$

For the problem (1) when we halve the step, the error get divided by 4, right?
That sounds reasonable, but how did you get it?
Seems to me you have to write it as a difference equation and solve it before you can say anything...

#### mathmari

##### Well-known member
MHB Site Helper
That sounds reasonable, but how did you get it?
Seems to me you have to write it as a difference equation and solve it before you can say anything...
The solution of the original differential is given. It's $y(x)=sin(\pi x)$.
Also the functions $f,q$ are given, $f(x)=\pi^4 \sin(\pi x)-\pi^2 q(x) \sin(\pi x)$, $q=\sin(x)^2$.

The finite difference method for $w$ is:
$$\frac{W_{j-1}-2W_j+W_{j+1}}{h^2}-q(x_j)W_j=f(x_j), j=1,...,J-1$$
$$W_0=W_J=0$$

where $W_j$ is the approximation of $w(x_j)$.

The finite difference method for $y$ is:
$$\frac{Y_{j-1}-2Y_j+Y_{j+1}}{h^2}=W_j, j=1,...,J-1$$
$$Y_0=Y_J=0$$

where $Y_j$ is the approximation of $y(x_j)$.

In an other exercise there was a finite difference method for a second order problem of two points, and when I halved the step, the error got divided by 4. So I thought that so it would be for all second order problems.

#### Klaas van Aarsen

##### MHB Seeker
Staff member
The solution of the original differential is given. It's $y(x)=sin(\pi x)$.
Also the functions $f,q$ are given, $f(x)=\pi^4 \sin(\pi x)-\pi^2 q(x) \sin(\pi x)$, $q=\sin(x)^2$.

The finite difference method for $w$ is:
$$\frac{W_{j-1}-2W_j+W_{j+1}}{h^2}-q(x_j)W_j=f(x_j), j=1,...,J-1$$
$$W_0=W_J=0$$

where $W_j$ is the approximation of $w(x_j)$.

The finite difference method for $y$ is:
$$\frac{Y_{j-1}-2Y_j+Y_{j+1}}{h^2}=W_j, j=1,...,J-1$$
$$Y_0=Y_J=0$$

where $Y_j$ is the approximation of $y(x_j)$.

In an other exercise there was a finite difference method for a second order problem of two points, and when I halved the step, the error got divided by 4. So I thought that so it would be for all second order problems.
Can it be that it should be $f(x)=\pi^4 \sin(\pi x)+\pi^2 q(x) \sin(\pi x)$?
With a $+$ instead of a $-$?

Assuming it is and rewriting the difference equation, I get:
$$W_{j+1}=h^2f(x_j) + (2 + h^2 q(x_j))W_j-W_{j-1}$$

If I solve it for h=0.1, h=0.05, and h=0.025, I get increasing errors.
I think it is numerically unstable.

That makes sense because an error in $W_j$ is blown up by at least a factor 2.

#### mathmari

##### Well-known member
MHB Site Helper
Can it be that it should be $f(x)=\pi^4 \sin(\pi x)+\pi^2 q(x) \sin(\pi x)$?
With a $+$ instead of a $-$?
No, it is like I wrote it, with a $-$.

Last edited:

#### Klaas van Aarsen

##### MHB Seeker
Staff member
No, it is like I wrote it, with a $-$.
Well, let's verify.
We have:
\begin{array}{}
y(x)&=&sin(\pi x) & (1)\\
f(x)&=&\pi^4 \sin(\pi x)-\pi^2 q(x) \sin(\pi x) & (2)\\
q(x)&=&\sin(x)^2 & (3)\\
w(x)&=&y''(x) & (4) \\
w''(x)-q(x)w(x)&=&f(x) & (5)
\end{array}

Substitute (1) in (4) to get:
$$w(x) = -\pi^2\sin(\pi x) \qquad\qquad\qquad (6)$$
Another 2 derivatives to get:
$$w''(x) = \pi^4\sin(\pi x) \qquad\qquad\qquad (7)$$
Substitute (3), (6) and (7) in (5) to get:
$$w''(x)-q(x)w(x) = \pi^4\sin(\pi x) - \sin(x)^2 \cdot -\pi^2\sin(\pi x) = \pi^4\sin(\pi x) + \pi^2 \sin(x)^2 \sin(\pi x)$$
That looks like $f(x)$ except for the $+$ sign...
It looks like either $f(x)$ is wrong or your differential equation is wrong...

#### mathmari

##### Well-known member
MHB Site Helper
Well, let's verify.
We have:
\begin{array}{}
y(x)&=&sin(\pi x) & (1)\\
f(x)&=&\pi^4 \sin(\pi x)-\pi^2 q(x) \sin(\pi x) & (2)\\
q(x)&=&\sin(x)^2 & (3)\\
w(x)&=&y''(x) & (4) \\
w''(x)-q(x)w(x)&=&f(x) & (5)
\end{array}

Substitute (1) in (4) to get:
$$w(x) = -\pi^2\sin(\pi x) \qquad\qquad\qquad (6)$$
Another 2 derivatives to get:
$$w''(x) = \pi^4\sin(\pi x) \qquad\qquad\qquad (7)$$
Substitute (3), (6) and (7) in (5) to get:
$$w''(x)-q(x)w(x) = \pi^4\sin(\pi x) - \sin(x)^2 \cdot -\pi^2\sin(\pi x) = \pi^4\sin(\pi x) + \pi^2 \sin(x)^2 \sin(\pi x)$$
That looks like $f(x)$ except for the $+$ sign...
It looks like either $f(x)$ is wrong or your differential equation is wrong...
I asked the professor and he told me that the function $f$ was wrong. It is with a $+$.

#### Klaas van Aarsen

##### MHB Seeker
Staff member
I asked the professor and he told me that the function $f$ was wrong. It is with a $+$.
Good!
So do you know how to solve the difference equation?

#### mathmari

##### Well-known member
MHB Site Helper
Good!
So do you know how to solve the difference equation?
I implemented the method to find the approximations $Y_j$ and found the following errors:
for $J=5$ $\text{error}=0.064097$
for $J=10$ $\text{error}=0.016395$

Could you tell me if these errors are right?

#### Klaas van Aarsen

##### MHB Seeker
Staff member
I implemented the method to find the approximations $Y_j$ and found the following errors:
for $J=5$ $\text{error}=0.064097$
for $J=10$ $\text{error}=0.016395$

Could you tell me if these errors are right?
Which method did you implement?
What is the reference point that you are using for your errors?

What I did, is use the shooting method to find the solution to your difference equation.
I compared that with the original function.
My error is the largest difference between the two (that occurs in the middle of the interval).
Chances are that you effectively did the same...

#### mathmari

##### Well-known member
MHB Site Helper
Which method did you implement?
The finite difference method.

What is the reference point that you are using for your errors?

What I did, is use the shooting method to find the solution to your difference equation.
I compared that with the original function.
My error is the largest difference between the two (that occurs in the middle of the interval).
Chances are that you effectively did the same...
I wrote a program in C, that calculates first the approximations $W_j$, for each $j$, using the Gauss elimination for tridiagonal systems. Then having found these approximations, it calculates the approximations $Y_j$ with the same way... The error is $max_{1 \leq j \leq J-1}|Y_j-y(x_j)|$, where $y$ is the real solution.

#### Klaas van Aarsen

##### MHB Seeker
Staff member
The finite difference method.

I wrote a program in C, that calculates first the approximations $W_j$, for each $j$, using the Gauss elimination for tridiagonal systems. Then having found these approximations, it calculates the approximations $Y_j$ with the same way... The error is $max_{1 \leq j \leq J-1}|Y_j-y(x_j)|$, where $y$ is the real solution.
I have for J=10 an error in y of 0.1514.

#### mathmari

##### Well-known member
MHB Site Helper
I have for J=10 an error in y of 0.1514.
Oh.. Ok.. I will check my program... And for J=20? When you double J, what happens with the error?

#### Klaas van Aarsen

##### MHB Seeker
Staff member
Oh.. Ok.. I will check my program... And for J=20? When you double J, what happens with the error?
With J=20, the maximum error in y that I get is 0.67848.

#### mathmari

##### Well-known member
MHB Site Helper
With J=20, the maximum error in y that I get is 0.67848.
Maybe I've found wrong the approximations..

For $J=4$ the approximations $W_j$ are:
69.305796
-10.379758
-7.337134

and the approximations $Y_j$:
69.305796
-1.2877786
-0.414607

The first one of both approximations isn't right, is it?
Could tell me if your result are different?

#### Klaas van Aarsen

##### MHB Seeker
Staff member
Maybe I've found wrong the approximations..

For $J=4$ the approximations $W_j$ are:
69.305796
-10.379758
-7.337134

and the approximations $Y_j$:
69.305796
-1.2877786
-0.414607

The first one of both approximations isn't right, is it?
Could tell me if your result are different?
My previous result for $J=20$ was wrong, the error in $Y$ should have been $0.05877$.

With $J=4$, I'm getting for $W$:
-6.708713243
-7.283950347
-3.563132688

And for $Y$:
0.452783069
0.450319241
0.225159621

for an error in $Y$ of $0.54968$.

#### mathmari

##### Well-known member
MHB Site Helper
My previous result for $J=20$ was wrong, the error in $Y$ should have been $0.05877$.

With $J=4$, I'm getting for $W$:
-6.708713243
-7.283950347
-3.563132688

And for $Y$:
0.452783069
0.450319241
0.225159621

for an error in $Y$ of $0.54968$.
I cannot find my error..

To find the approximations $W_j$ I did the following:

Code:
for(j=0; j<J-1; j++){
x=(j+1)*h;
b[j]=f(x);
}
For the Gauss elimination for tridiagonal systems, the system is $Aw=b$, where $b$ is calculated above and the matrix $A$ is tridiagonal which has the form:
Diagonal:-2.0/(h*h)-q(x), where x=(j+1)*h
Upper diagonal:1.0/(h*h)
Lower diagonal:1.0/(h*h)

Do you see something wrong? Did you find these approximations with an other way?

#### Klaas van Aarsen

##### MHB Seeker
Staff member
I cannot find my error..

To find the approximations $W_j$ I did the following:

Code:
for(j=0; j<J-1; j++){
x=(j+1)*h;
b[j]=f(x);
}
For the Gauss elimination for tridiagonal systems, the system is $Aw=b$, where $b$ is calculated above and the matrix $A$ is tridiagonal which has the form:
Diagonal:-2.0/(h*h)-q(x), where x=(j+1)*h
Upper diagonal:1.0/(h*h)
Lower diagonal:1.0/(h*h)

Do you see something wrong? Did you find these approximations with an other way?
That looks correct.
However, your value of 69.305796 for $W_1$ is clearly way too large.

I used a different method: the shooting method.
I wrote each value as a recursive relation depending on the previous 2 values.
Then I changed $W_1$ just as long until $W_J$ became 0.

Either way, you should get the same values.
How did you solve $Aw=b$?
Did you verify by evaluating $Aw$ with the result for $w$ you found?

#### mathmari

##### Well-known member
MHB Site Helper
However, your value of 69.305796 for $W_1$ is clearly way too large.
Yes...I ran the program for several $J$, and for $J>4$ the values for $W_j$ are not so large...
For example for $J=5$ I get:
for $W_j$:
-5.671676
-8.580075
-7.553415
-3.738234

for $Y_j$:
0.538176
0.849485
0.817591
0.483560

How did you solve $Aw=b$?
A is a matrix of dimension (J-1)x3.
It has the form:
$\begin{bmatrix} 0 & d_1 & u_1\\ l_1 & d_2 &u_2 \\ l_2 & d_3 &u_3 \\ ... & ... & ...\\ l_{J_3} & d_{J-2} & u_{J-2}\\ l_{J-2} & d_{J-1} & 0 \end{bmatrix}$
where $l_1,...,l_{J-1}$ are the elements of the lowdiagonal, $u_1,...,u_{J-1}$ the elements of the upperdiagonal and $d_1,...,d_{J-1}$ the elements of the diagonal.
then we shift the elements of the first row one place to the left and then we implement the Gauss elimination with pivot. We do that at each step..

Did you verify by evaluating $Aw$ with the result for $w$ you found?
No, I haven't done it yet..

#### Klaas van Aarsen

##### MHB Seeker
Staff member
A is a matrix of dimension (J-1)x3.
It has the form:
$\begin{bmatrix} 0 & d_1 & u_1\\ l_1 & d_2 &u_2 \\ l_2 & d_3 &u_3 \\ ... & ... & ...\\ l_{J_3} & d_{J-2} & u_{J-2}\\ l_{J-2} & d_{J-1} & 0 \end{bmatrix}$
where $l_1,...,l_{J-1}$ are the elements of the lowdiagonal, $u_1,...,u_{J-1}$ the elements of the upperdiagonal and $d_1,...,d_{J-1}$ the elements of the diagonal.
then we shift the elements of the first row one place to the left and then we implement the Gauss elimination with pivot. We do that at each step..
That sounds about right... but it's easy to make a mistake there.

No, I haven't done it yet..
So will you?

#### mathmari

##### Well-known member
MHB Site Helper
So will you?
For $J=5$ the matrix $A$ is $\begin{bmatrix} -50.039470& 25 & 0 & 0\\ 25& -50.151647 & 25 &0 \\ 0 & 25 &-50.318821 & 25\\ 0& 0 & 25& -50.514600 \end{bmatrix}$

and $b=\begin{bmatrix} 57.484598\\ 94.064990 \\ 95.634182\\ 60.240927 \end{bmatrix}$

Using the Matlab command w=A\b, I get the same $w$ as I get with my program:
$w=\begin{bmatrix} -5.992187\\ -9.694450 \\ -9.692919\\ -5.989633 \end{bmatrix}$

#### mathmari

##### Well-known member
MHB Site Helper
That means that the solution of the system is right, isn't it ?