How do I prepare the following ODE for ode45?
1 次查看(过去 30 天)
显示 更早的评论
Sergio Manzetti
2018-2-20
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 个评论
Sergio Manzetti
2018-2-20
编辑:Sergio Manzetti
2018-2-20
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?
采纳的回答
Torsten
2018-2-20
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 个评论
Sergio Manzetti
2018-2-21
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
2018-2-21
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
2018-2-22
编辑:Sergio Manzetti
2018-2-22
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
2018-2-22
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
2018-2-22
Thanks Torsten.
So ODE45 represents the form y''' + y(1-x^2)=0 as
[ y(2) y(3) y(1)*(2+x^2)]
Torsten
2018-2-22
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
2018-2-22
编辑:Sergio Manzetti
2018-2-22
Thanks! I did it read it, but haven't done this in 2 years, so it was not clear.
Cheers
Sergio Manzetti
2018-2-22
编辑:Sergio Manzetti
2018-2-22
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
2018-2-22
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
2018-2-22
编辑:Sergio Manzetti
2018-2-22
Thanks Torsten! PS: Is it possible with this to get the third order derivative of the solution?
Jan
2018-2-22
@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
2018-2-28
Hi Jan, is it possible to plot the square modulus of the numerical solution following your suggestion here ?
更多回答(1 个)
Sergio Manzetti
2018-2-28
编辑:Sergio Manzetti
2018-2-28
(abs(y))^2
if y is the solution
11 个评论
Sergio Manzetti
2018-3-9
编辑:Sergio Manzetti
2018-3-9
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
2018-3-9
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
2018-3-9
编辑:Sergio Manzetti
2018-3-9
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
2018-3-9
编辑:Sergio Manzetti
2018-3-9
OK, this is quite a new way to think...so one lists the levels of derivation of y as such?
Torsten
2018-3-9
You transform a higher order ODE to a system of first-order ODEs.
I already gave you the link to digest this.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
标签
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)