Solving the Initial Value Problem for a Wave Using the Forward Euler Method

[evinda]

Hello! (Wave)

We consider the initial value problem

$$x'(t)=-y(t), t \in [0,1] \\ y'(t)=x(t), t \in [0,1] \\ x(0)=1, y(0)=0$$

I want to solve approximately the above problem using the forward Euler method in uniform partition of 100 and 200 points.

I have written the following code in matlab:

N=100;
h=1/N;
y=zeros(N);
A=[0 -1;1 0];
for (i=1:100)
    y=(eye(N,N)+A*h)*y
    i=i+1;
end
y

Is my code right? (Thinking)

 

[I like Serena]

Hey evinda! (Smile)

The dimensions of A do not match eye(N,N), nor y.
Shouldn't it be eye(2)?

And where is the initial condition $x(0)=1, y(0)=0$?
Shouldn't we have y=[1;0] instead of y=zeroes(N)? (Wondering)

 

[HOI]

You can check your work by solving the equations analytically. The equations are
x′(t)=−y(t) and y′(t)=x(t). Differentiating the first equation again, x'(t)= -y'(t)= -x. That is a linear ordinary differential equation with constant coefficients. The general solution is x(t)= A cos(t)+ B sin(t). y(t)= -x'(t)= A sin(t)- B cos(t).

The initial condition is x(0)= A= 1, y(0)= -B= 0.

x(t)= cos(t), y(t)= sin(t).
Check your values against that.

 

[evinda]

Hey evinda! (Smile)

The dimensions of A do not match eye(N,N), nor y.
Shouldn't it be eye(2)?

And where is the initial condition $x(0)=1, y(0)=0$?
Shouldn't we have y=[1;0] instead of y=zeroes(N)? (Wondering)

Ah I see... So the correct code would be the following, right? (Blush)

N=1/100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
for (i=1:N)
      y=(eye(2)+A*h)*y;
      i=i+1;
end
y

Using the trapezoid method, the code would get the following form, right?

N=100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
for (i=1:N)
     y=inv(eye(2)-h/2*A)*(eye(2)+h/2*A)*y;
     i=i+1;
end
y

 

[I like Serena]

Ah I see... So the correct code would be the following, right? (Blush)

N=1/100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
for (i=1:N)
      y=(eye(2)+A*h)*y;
      i=i+1;
end
y

The [M]i=i+1[/M] is redundant, but yes, that will give the Nth approximation of the Forward Euler method. (Nod)

Using the trapezoid method, the code would get the following form, right?

N=100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
for (i=1:N)
     y=inv(eye(2)-h/2*A)*(eye(2)+h/2*A)*y;
     i=i+1;
end
y

How did you get that? (Wondering)

If I apply Heun's method, also known as the Explicit trapezoidal rule, I get:
$$y_{i+1} = y+\frac h2\big[Ay+A(y+hAy)\big] = y+hAy + \frac{h^2}{2}A^2y$$

 

[evinda]

The [M]i=i+1[/M] is redundant, but yes, that will give the Nth approximation of the Forward Euler method. (Nod)

Ah yes, nice! (Smirk)

If I apply Heun's method, also known as the Explicit trapezoidal rule, I get:
$$y_{i+1} = y+\frac h2\big[Ay+A(y+hAy)\big] = y+hAy + \frac{h^2}{2}A^2y$$

Ok. Then the code should look as follows, shouldn't it? (Thinking)

N=100;
h=1/N;
y=[1;0];
A=[0 -1; 1 0];
for (i=1:N)
      y=(eye(2)+h*A+h^2/2*A^2)*y;
end
y

 

[I like Serena]

Yep.
And it's already pretty close to [M][cos(1);sin(1)][/M]. (Happy)

 

[evinda]

Yep.
And it's already pretty close to [M][cos(1);sin(1)][/M]. (Happy)

You mean the solution y? How do we deduce this?

Using the Runge-Kutta method of fourth order, we have the following code, right?

N=100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
for (i=1:N)
      k1=A*y
      k2=A*(y+1/2*k1*h);
      k3=A*(y+1/2*k2*h);
      k4=A*(y+k3*h);
      y=y+1/6*(k1+2*k2+2*k3+k4)*h;
end
y

 

[I like Serena]

You mean the solution y? How do we deduce this?

It's the solution to the initial value problem, isn't it?
It we substitute $x(t)=\cos t$ and $y(t)=\sin t$, it satisfies the conditions. (Thinking)

Using the Runge-Kutta method of fourth order, we have the following code, right?

N=100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
for (i=1:N)
      k1=A*y
      k2=A*(y+1/2*k1*h);
      k3=A*(y+1/2*k2*h);
      k4=A*(y+k3*h);
      y=y+1/6*(k1+2*k2+2*k3+k4)*h;
end
y

Yep. (Nod)

 

[evinda]

It's the solution to the initial value problem, isn't it?
It we substitute $x(t)=\cos t$ and $y(t)=\sin t$, it satisfies the conditions. (Thinking)

Yes, it is. (Nod)

Yep. (Nod)

Nice... (Happy)

I also want to solve the problem using the Gauss-Legendre method, $q=2$, in uniform partition of 100 and 200 points.

Which is the formula of this method? I have looked for it in the internet, but I didn't find something helpful. (Nerd)

 

[I like Serena]

I also want to solve the problem using the Gauss-Legendre method, $q=2$, in uniform partition of 100 and 200 points.

Which is the formula of this method? I have looked for it in the internet, but I didn't find something helpful. (Nerd)

Until now, I only knew of the Gauss–Legendre quadrature for integration.
Anyway, http://www.wseas.us/e-library/conferences/2006bucharest/papers/518-127.pdf?origin=publication_detail seems to explain how to use it for IVP's as a variant of Runge-Kutta, including a dedicated 2-point formula. (Thinking)

 

[evinda]

Until now, I only knew of the Gauss–Legendre quadrature for integration.
Anyway, http://www.wseas.us/e-library/conferences/2006bucharest/papers/518-127.pdf?origin=publication_detail seems to explain how to use it for IVP's as a variant of Runge-Kutta, including a dedicated 2-point formula. (Thinking)

So the code will have the following form, right? (Thinking)

N=100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
a1=0.5*(1-1/sqrt(3));
a2=0.5*(1+1/sqrt(3));
for (i=1:N)
      k1=(A+a1*h*A^2)*y;
      k2=A*(y+a2*h*k1);
      y=y+h/2*(k1+k2);
end
y

 

[I like Serena]

So the code will have the following form, right? (Thinking)

N=100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
a1=0.5*(1-1/sqrt(3));
a2=0.5*(1+1/sqrt(3));
for (i=1:N)
      k1=(A+a1*h*A^2)*y;
      k2=A*(y+a2*h*k1);
      y=y+h/2*(k1+k2);
end
y

Yep. (Nod)

 

[evinda]

Yep. (Nod)

Great. (Blush) I also want to represent in common graphs the corresponding values of $(x^n)^2+(y^n)^2$.
So do we add at each of the above for-loops the command [M] plot(y(1)^2+y(2)^2)[/M] after the calculation of y? (Thinking)

 

[I like Serena]

Great. (Blush) I also want to represent in common graphs the corresponding values of $(x^n)^2+(y^n)^2$.
So do we add at each of the above for-loops the command [M] plot(y(1)^2+y(2)^2)[/M] after the calculation of y? (Thinking)

Well, according to description of plot, we have:

plot(Y) creates a 2-D line plot of the data in Y versus the index of each value.​

That means that the first time we get a plot of the first point at x-coordinate 1.
The second time that plot is replaced by a plot of a single point that is again at x-coordinate 1.
And so on, ending in the plot of a single point that corresponds to the last iteration.
That's not what we want, is it? (Worried)

To add to the plot instead of replacing it, we can call [M]hold on[/M].
And perhaps we should specify an x-coordinate of, say, [M]i*h[/M]? (Wondering)

 

[evinda]

To add to the plot instead of replacing it, we can call [M]hold on[/M].
And perhaps we should specify an x-coordinate of, say, [M]i*h[/M]? (Wondering)

I haven't really understood how we use the command plot. Could you maybe explain it further to me? (Worried)

 

[I like Serena]

I haven't really understood how we use the command plot. Could you maybe explain it further to me? (Worried)

Suppose we do:

N=100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
plot(0, y(1)^2+y(2)^2);
hold on;
for (i=1:N)
      y=(eye(2)+A*h)*y;
      plot(i*h, y(1)^2+y(2)^2);
end
hold off;
y

We use the plot(X,Y) command, which is described as:

plot(X,Y) creates a 2-D line plot of the data in Y versus the corresponding values in X.​

Then the first plot command will create the plot and put one value in there at x-coordinate 0.
Calling [M]hold on[/M] means that we want to hold the plot, such that subsequent plot commands get added to the same plot until we call [M]hold off[/M]. Without it, the previous plot command would be discarded.
And [M]plot(i*h, y(1)^2+y(2)^2)[/M] plots the desired value at the corresponding x-coordinate, which is [M]x=i*h[/M].

The result is:
View attachment 8104
(Thinking)

 

[I like Serena]

Alternatively, we could first gather all the data, and afterwards create the plot.
For instance with:

N=100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
tt=[0];
yy=[y];
for (i=1:N)
      y=(eye(2)+A*h)*y;
      tt=[tt, i*h];
      yy=[yy, y];
end
plot(tt, yy(1,:).^2+yy(2,:).^2);

The advantage is that we can easily add other graphs in the same plot.
For instance with the actual x and y values.
Or with the graphs for the other approximations. (Thinking)

 

[evinda]

Suppose we do:

N=100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
plot(0, y(1)^2+y(2)^2);
hold on;
for (i=1:N)
      y=(eye(2)+A*h)*y;
      plot(i*h, y(1)^2+y(2)^2);
end
hold off;
y

We use the plot(X,Y) command, which is described as:

plot(X,Y) creates a 2-D line plot of the data in Y versus the corresponding values in X.​

Then the first plot command will create the plot and put one value in there at x-coordinate 0.
Calling [M]hold on[/M] means that we want to hold the plot, such that subsequent plot commands get added to the same plot until we call [M]hold off[/M]. Without it, the previous plot command would be discarded.
And [M]plot(i*h, y(1)^2+y(2)^2)[/M] plots the desired value at the corresponding x-coordinate, which is [M]x=i*h[/M].

The result is:

(Thinking)

Ah I see... Like that we get the plot of the approximation of $x^2+y^2$ at each time $t_i=0+ih$, right?

- - - Updated - - -

Alternatively, we could first gather all the data, and afterwards create the plot.
For instance with:

N=100;
h=1/N;
y=[1;0];
A=[0 -1;1 0];
tt=[0];
yy=[y];
for (i=1:N)
      y=(eye(2)+A*h)*y;
      tt=[tt, i*h];
      yy=[yy, y];
end
plot(tt, yy(1,:).^2+yy(2,:).^2);

The advantage is that we can easily add other graphs in the same plot.
For instance with the actual x and y values.
Or with the graphs for the other approximations. (Thinking)

Could you explain to me the commands

[M]tt=[0];
yy=[y];[/M]

and

[M]tt=[tt, i*h];
yy=[yy, y]; [/M]

?

With tt=[0]; for example, how do we know that the 2-dimensionsal zero vector is meant?

 

[I like Serena]

Ah I see... Like that we get the plot of the approximation of $x^2+y^2$ at each time $t_i=0+ih$, right?

Indeed. (Nod)

Could you explain to me the commands

[M]tt=[0];
yy=[y];[/M]

and

[M]tt=[tt, i*h];
yy=[yy, y]; [/M]

?

With tt=[0]; for example, how do we know that the 2-dimensionsal zero vector is meant?

It's not. (Shake)

[M]tt=[0][/M] creates a row vector of t-values with the single scalar value 0 in it.
[M]yy=[y][/M] creates a matrix of y-vectors with with the single column vector y in it.
[M]tt=[tt, i*h][/M] creates a new row vector by adding [M]i*h[/M] to the right of [M]tt[/M].
[M]yy=[yy, y][/M] creates a new matrix by appending the column vector [M]y[/M] to the right of the matrix [M]yy[/M].

Afterwards [M]yy[/M] contains all the iterations of [M]y[/M], and [M]tt[/M] contains the corresponding t-values. (Thinking)

 

[evinda]

Indeed. (Nod)
It's not. (Shake)

[M]tt=[0][/M] creates a row vector of t-values with the single scalar value 0 in it.
[M]yy=[y][/M] creates a matrix of y-vectors with with the single column vector y in it.
[M]tt=[tt, i*h][/M] creates a new row vector by adding [M]i*h[/M] to the right of [M]tt[/M].
[M]yy=[yy, y][/M] creates a new matrix by appending the column vector [M]y[/M] to the right of the matrix [M]yy[/M].

Afterwards [M]yy[/M] contains all the iterations of [M]y[/M], and [M]tt[/M] contains the corresponding t-values. (Thinking)

Nice... Thank you very much! (Happy)