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.
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.
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:
$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?
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?
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?
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 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}$
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$?
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 $$\frac{1}{Dt^2}+\frac{2}{h^2}$$ and for their neighbors $$- \frac{1}{h^2}$$.
Hmm, you must have rearranged it differently from how I did it.
I rearranged the formula you gave earlier and found for the diagonal elements $$\frac{1}{Dt^2}+\frac{2}{h^2}$$ and for their neighbors $$- \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?
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)$
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)$
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] ?
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.
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?
I read again the exercise and I saw that we need an array $A$ with $J-1$ elements since the zero elements don't have to be saved. That means that I have to define this array so: double A[J-1]; , or not?
So I changed the function that calculates $b$:
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-1; 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++;
}
C[J-2]=A[J-3]/(2*h*h)-A[J-2]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+2*B[J-2]/(Dt*Dt)+f(xx(J-1,J),tn(n,N));
pointer=&C[J-2];
printf("%lf\n",*pointer);
}
I read again the exercise and I saw that we need an array $A$ with $J-1$ elements since the zero elements don't have to be saved. That means that I have to define this array so: double A[J-1]; , or not?
So I changed the function that calculates $b$:
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-1; 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++;
}
C[J-2]=A[J-3]/(2*h*h)-A[J-2]*(h*h+Dt*Dt)/(Dt*Dt*h*h)+2*B[J-2]/(Dt*Dt)+f(xx(J-1,J),tn(n,N));
pointer=&C[J-2];
printf("%lf\n",*pointer);
}
Now there is no $A[J-1]$..Is it correct now?
Yep. That is correct.
Btw, any reason that the return value is "void *" instead of "void"?
And that you execute "pointer++" when you make no use of the incremented value?
