**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# How do I prepare the following ODE for ode45?

18 views (last 30 days)

Show older comments

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!

##### 2 Comments

Sergio Manzetti
on 20 Feb 2018

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?

### Accepted Answer

Torsten
on 20 Feb 2018

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.

##### 17 Comments

Sergio Manzetti
on 21 Feb 2018

Steven Lord
on 21 Feb 2018

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
on 22 Feb 2018

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?

Sergio Manzetti
on 22 Feb 2018

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 ?

Sergio Manzetti
on 22 Feb 2018

Thanks Torsten.

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

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

Torsten
on 22 Feb 2018

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
on 22 Feb 2018

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

Cheers

Sergio Manzetti
on 22 Feb 2018

Sergio Manzetti
on 22 Feb 2018

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

Sergio Manzetti
on 28 Feb 2018

### More Answers (1)

Sergio Manzetti
on 28 Feb 2018

Edited: Sergio Manzetti
on 28 Feb 2018

(abs(y))^2

if y is the solution

##### 11 Comments

Sergio Manzetti
on 9 Mar 2018

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!

Sergio Manzetti
on 9 Mar 2018

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?

Sergio Manzetti
on 9 Mar 2018

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

### See Also

### Categories

### Tags

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.