How do I prepare the following ODE for ode45?

Hello, I would like to solve the following ODE in ode45, but the example's on the site are not describing using higher order derivatives with non-linear terms.

The ODE is:

y''' = y(2+x^2)

initial conditions are: y(0)=0 y'(0)=0 y''(0)=0

Thanks!

[Torsten]

Google is your friend:

matlab & higher order odes

Best wishes

Torsten.

[Sergio Manzetti]

I got so far:

dYdX = @(X,Y) [Y(3) +  (x^2+2)*Y(1)];  % Differential equation

res = @(ya,yb) [ya(1);  ya(2);  yb(2)-1];                   % Boundary conditions

SolYinit = bvpinit([0 1E+1], [1; 1; 1]);

Fsol = bvp4c(dYdX, res, SolYinit);

X = Fsol.x;
F = Fsol.y;

figure(1)
plot(X, F)
legend('F_1', 'F_2', 'F_3', 3) 
grid

But the first line is not correct. Can you see what is missing?

Answer 1

[Torsten]

fun = @(x,y)[y(2);y(3);y(1)*(2+x^2)];
y0 = [0 0 0];
xspan = [0 5];
[X,Y] = ode45(fun,xspan,y0);
plot(X,Y(:,1),X,Y(:,2),X,Y(:,3));

Best wishes

Torsten.

[Sergio Manzetti]

Fantastic, thanks Torsten!

[Sergio Manzetti]

Hi Torsten, I am trying to expand the interval from 0 5 to -500 to +500, but I only get a flat line. Even when I go back to 0 5 in the interval, the graph again appears as a flat line...

[Steven Lord]

Your initial conditions state that y, its first derivative, and its second derivative are all 0 at x = 0. If you think about this from a physics standpoint where y is the position of your solution (let's think of this as a race car), y' is the velocity, and y'' is the acceleration this corresponds to you being at the starting line at a dead stop with your foot off the gas pedal.

Speaking very roughly, at the first time step (I'm treating x as time here) y'' is 0*something, so you never accelerate. If you don't accelerate you can't go faster than 0. If you don't go faster than 0, you can't move away from 0. Torsten's solution (correctly) demonstrates that the solution to your system of equations with that set of initial conditions is y(x) = 0.

Now if you were to tweak the initial conditions, say by stomping on the gas pedal as the race starts:

>> y0 = [0 0 1];
>> [X,Y] = ode45(fun,xspan,y0);
>> plot(X,Y(:,1),X,Y(:,2),X,Y(:,3));

You get something that looks a bit more interesting. For most of the time span your curves looks like they are at 0, but that's because of the scaling on the Y axis. If you zoom in you can see that y, its first derivative, and its second derivative take off very rapidly. [Given that the acceleration increases as the square of x, you get going REALLY quickly.]

[Sergio Manzetti]

Hi, unfortunately I have no information on what the acceleration is, but it can be given an arbitrary value to begin with. Thanks for this suggestion, I will work on this model instead. The first line was originally foreign to me, does it mean that y(0) = 0 , y'(0)=0 and y''(0)=1 ?

The plot becomes more interesting as the second derivative is set to 1, as you say. But how do I interpret the three lines here?

[Torsten]

y, y' and y'' (in this order).

Best wishes

Torsten.

[Sergio Manzetti]

Thanks, if I add a constant to the y''' term, it appears as:

fun = @(x,y)[y(2);y(3)*h;y(1)*(2+x^2)];
y0 = [0 0 1];
xspan = [0 5];
[X,Y] = ode45(fun,xspan,y0);
plot(X,Y(:,1),X,Y(:,2),X,Y(:,3));

where h is defined before. But why this weird structure in ODE45 to represent y''': y(2);y(3)*h ?

[Torsten]

It's not a "weird structure", but the usual way to reduce higher-order ODEs to a system of first-order ODEs. See e.g.

http://www.math.poly.edu/courses/ma2132/Notes/MA2132EquationsToSystems.pdf

If you define

fun = @(x,y)[y(2);y(3)*h;y(1)*(2+x^2)];

you solve e.g. y'''=h*y*(2+x^2).

Best wishes

Torsten.

[Sergio Manzetti]

Thanks Torsten.

So ODE45 represents the form y''' + y(1-x^2)=0 as

[ y(2) y(3) y(1)*(2+x^2)]

[Torsten]

No.

fun = @(x,y)[y(2);y(3);y(1)*(2+x^2)];

represents the system of ODEs

y1' = y2
y2' = y3
y3' = y1*(1-x^2)

And having a solution for y1, y2 and y3 of this system gives you a solution of y'''=y*(1-x^2) because y1'''=y2''=y3'=y1*(1-x^2).

Please read the document I linked to.

Best wishes

Torsten.

[Sergio Manzetti]

Thanks! I did it read it, but haven't done this in 2 years, so it was not clear.

Cheers

[Sergio Manzetti]

About the picture, I think I misunderstood you before, you said y, y' and y''. The colors are associated by default to y , y' and y''?

[Torsten]

The first curve is the solution of your third-order ODE.

Use

plot(X,Y(:,1))

if you only want to see this curve and not also its first and second derivative.

Best wishes

Torsten.

[Sergio Manzetti]

Thanks Torsten! PS: Is it possible with this to get the third order derivative of the solution?

[Jan]

@Sergio: If you want the 3rd derivative, simply use the already provided function:

fun = @(x,y)[y(2);y(3);y(1)*(2+x^2)];

After

[X,Y] = ode45(fun,xspan,y0);

you have X and Y (which is [y, y', y'']. So you can calculate the 3rd derivative easily:

d3y = y(:,1) .* (2 + X .^ 2);

[Sergio Manzetti]

Thanks Jan! Will give it a try!

[Sergio Manzetti]

Hi Jan, is it possible to plot the square modulus of the numerical solution following your suggestion here ?

[Torsten]

Maybe, if you tell us what the "square modulus of the numerical solution" is.

Answer 2

[Sergio Manzetti]

(abs(y))^2

if y is the solution

[Torsten]

SQM = abs(Y(:,1).^2)

[Sergio Manzetti]

Fantastic, Thanks Torsten!

[Torsten]

For Y real-valued,

SQM=Y(:,1).^2

suffices, of course.

[Sergio Manzetti]

Hi Torsten you wrote in the initial ODE45 command:

fun = @(x,y)[y(2);y(3);y(1)*(2+x^2)];
y0 = [1 0 0];
xspan = [0 5];
[X,Y] = ode45(fun,xspan,y0);
plot(X,Y(:,1),X,Y(:,2),X,Y(:,3));

but the first line fun appears to me as D3Y*(2+x^2)Y=0

Is this correct? If so it is not the function I wrote.

Sorry I have to re-check.

What I wanted to calculate is

y''' = y(2+x^2)

thus D3Y - y(2+x"2)=0

so would this code reflect that?

fun = @(x,y)[y(2);y(3);y(1)-(2+x^2)];
  y0 = [1 0 0];
  xspan = [0 5];
  [X,Y] = ode45(fun,xspan,y0);
  plot(X,Y(:,1),X,Y(:,2),X,Y(:,3));

Thanks!

[Torsten]

The first line is

(1) y1'(x) = y2(x)
(2) y2'(x) = y3(x)
(3) y3'(x) = y1(x)*(2+x^2).

Thus the equation for y1 is

y1(x)*(2+x^2) = (from (3)) y3'(x) = (from (2)) y2''(x) = (from (1)) y1'''(x)

which means that y1 is the function you searched for.

[Sergio Manzetti]

Thanks Torsten, does this mean that the third derivative in ODE45 is given by:

(1) y1'(x) = y2(x)
(2) y2'(x) = y3(x)

so y(2);y(3)?

Should I want to multiply a constant to the D3Y part, is it multiplied as such?:

fun = @(x,y)[y(2);y(3) *h;y(1)*(2+x^2)];

where h is the constant?

[Torsten]

y2 = y'
y3 = y''

[Sergio Manzetti]

where is y''' ?

y''' = y(2+x^2)

[Torsten]

y''' = y3' = y1*(2+x^2)

That's (3).

[Sergio Manzetti]

OK, this is quite a new way to think...so one lists the levels of derivation of y as such?

[Torsten]

You transform a higher order ODE to a system of first-order ODEs.

I already gave you the link to digest this.