MHB Finite difference method-convergence

[mathmari]

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}?

 

[mathmari]

I would also know how the error changes while the step size is getting smaller? Is it getting smaller or bigger??

 

[I like Serena]

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}}$.
Can you provide some more information?

 

[mathmari]

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??

 

[I like Serena]

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]

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.

 

[I like Serena]

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]

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 $-$.

 

[I like Serena]

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, 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 $+$.

 

[I like Serena]

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]

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?

 

[I like Serena]

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]

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.

 

[I like Serena]

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]

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?

 

[I like Serena]

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]

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?

 

[I like Serena]

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]

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:

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?

 

[I like Serena]

I cannot find my error.. :(

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

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]

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..

 

[I like Serena]

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? (Wasntme)

 

[mathmari]

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]

That means that the solution of the system is right, isn't it ?

 

[I like Serena]

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 $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}$

That means that the solution of the system is right, isn't it ?

Yep.
Did you fix something?
Because in your previous post you had a different solution for w...

 

[mathmari]

Yep.
Did you fix something?
Because in your previous post you had a different solution for w...

At the following part I replaced $J-1$ with $J$.

for(j=0; j<J-1; j++){
                x=(j+1)*h;
		b[j]=f(x);	
}

for(j=0; j<J; j++){
                x=(j+1)*h;
		b[j]=f(x);	
}

 

[I like Serena]

At the following part I replaced $J-1$ with $J$.

for(j=0; j<J; j++){
                x=(j+1)*h;
		b[j]=f(x);	
}

Oh?
That should not make a difference??
Can there be a $\pm 1$ mistake afterwards?

 

[mathmari]

Can there be a $\pm 1$ mistake afterwards?

What do you mean?

 

[I like Serena]

What do you mean?

You have measurements at $x_0, x_1, ..., x_J$.
That are $(J+1)$ measurements.
However, you can leave out $x_0$ and $x_J$ since they are both given as boundary conditions: they are 0.
The leaves you with $(J-1)$ measurements.
So your matrix should be $(J-1) \times (J-1)$, your $b$ should have $(J-1)$ entries, and your resulting $w$ should have $(J-1)$ entries.

However, your "fix" gives you $J$ entries instead of $(J-1)$ entries.
That is $1$ off.
If it helps you to get the right answer, it means that later on you are also $1$ off.
So before you were probably calculating with an uninitialized variable.

 

[mathmari]

You have measurements at $x_0, x_1, ..., x_J$.
That are $(J+1)$ measurements.
However, you can leave out $x_0$ and $x_J$ since they are both given as boundary conditions: they are 0.
The leaves you with $(J-1)$ measurements.
So your matrix should be $(J-1) \times (J-1)$, your $b$ should have $(J-1)$ entries, and your resulting $w$ should have $(J-1)$ entries.

However, your "fix" gives you $J$ entries instead of $(J-1)$ entries.
That is $1$ off.
If it helps you to get the right answer, it means that later on you are also $1$ off.
So before you were probably calculating with an uninitialized variable.

Yes,you are right!

I looked again my code and realized that I didn't fix that $J$.
I have written the following in a function:

fun(b,J){
         for(j=0; j<J-1; j++)
		b[j]=f(X(j,J));	
         
 }

At the beginning I called the function fun(b,J-1) and that is wrong, since X is calculated in the function

X(j,J){
      x=(j+1)/J;
 }

Then I fixed it by calling the function by fun(b,J).

 

[I like Serena]

Yes,you are right!

I looked again my code and realized that I didn't fix that $J$.
I have written the following in a function, fun(b,J):

for(j=0; j<J-1; j++){
                x=(j+1)*h;
		b[j]=f(x);	
}

At the beginning I called the function fun(b,J-1) and in the function it was for(j=0; j<J; j++) and that is wrong.
Then I fixed it by calling the function by fun(b,J).

Aha!