Welcome to our community

Be a part of something great, join today!

Finite difference method-convergence

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
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 [tex] \frac{log(\frac{e_{1}}{e_{2}})}{log(\frac{J_{2}}{J_{1}})} [/tex], right? where [tex] e_{i}=max|y_{j}^{J_{i}}-y(x_{j}^{J_{i}})| [/tex], [tex] 0<=j<=J_{i} [/tex]. How can I find the [tex] J_{1} [/tex] and [tex] J_{2} [/tex]?
 
Last edited:

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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 [tex] \frac{log(\frac{e_{1}}{e_{2}})}{log(\frac{J_{2}}{J_{1}})} [/tex], right? where [tex] e_{i}=max|y_{j}^{J_{i}}-y(x_{j}^{J_{i}})| [/tex], [tex] 0<=j<=J_{i} [/tex]. How can I find the [tex] J_{1} [/tex] and [tex] J_{2} [/tex]?
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}}$.
Can you provide some more information?
 

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
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}}$.
Can you provide some more information?
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
Mar 5, 2012
8,780
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
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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
Apr 14, 2013
4,047
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
Mar 5, 2012
8,780
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... :eek:
 

mathmari

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

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
Apr 14, 2013
4,047
Which method did you implement?
The finite difference method.

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

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
Mar 5, 2012
8,780
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
Apr 14, 2013
4,047

Klaas van Aarsen

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


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

mathmari

Well-known member
MHB Site Helper
Apr 14, 2013
4,047
So will you? (Wasntme)
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
Apr 14, 2013
4,047
That means that the solution of the system is right, isn't it ?