ODE45 Code Application

Question

James 2020-6-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/539383-ode45-code-application

评论： James 2020-6-1

M = 100/6*pi*0.04^3; g=9.8; A=0.25*pi*0.04^2;
Re = @(U) 1.25*0.04*U/0.000015;
CD = @(Re) 24./Re+(2.6*Re/5)./(1+(Re/5).^1.52)+(0.411*(Re/263000).^(-7.94))./(1+(Re/263000).^(-8))+(0.25*Re/10^6)./(1+Re/10^6);
function Udot = Velocity(t,U)
    Udot = g -(0.5*1.25*U^2*A*CD)/M;
    [t,U]=ode45 (@Velocity, [0, 20], 0);
    plot(t,U)
end

I'm in agony with this coding for 5hours..... I can solve easy function but it is very complicate function for me. I'm confused how to make up 'Re' 'CD' coding clearly.

Thanks for reading.. how could I get right plotting?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Goodmanson 2020-6-1

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/539383-ode45-code-application#answer_443963

编辑：David Goodmanson 2020-6-1

Hi James,

[t,U] = ode45(@Velocity, [0, 20], .001);
plot(t,U)
function Udot = Velocity(t,U)
M = 100/6*pi*0.04^3;
g=9.8;
A=0.25*pi*0.04^2;
Re = 1.25*0.04*U/0.000015;
CD  = 24./Re+(2.6*Re/5)./(1+(Re/5).^1.52)+(0.411*(Re/263000).^(-7.94))./(1+(Re/263000).^(-8))+(0.25*Re/10^6)./(1+Re/10^6);
Udot = g -(0.5*1.25*U^2*A*CD)/M;
end

You need to get ode45 outside of the functon definition. And you can calculate Re and CD without defining functions for them. If you do define Re and CD as functions, then calculating CD would require using Re(U) everwhere instead of just Re, and similarly CD would have to supplied with an argument. Easier here to not do that.

I changed the initial value of U to .001 since the 24/Re term blows up if you use zero. Changing the initial value does not seem to make much difference as long as it's small.