Ordinary differential equation - pratical problem

2016 7 9
0 个回答
更新时间：2016 7 13
7 次查看（30 天）

0 个投票

Hello,
I have a single equation of the form: u'=a*u^2+b*u^1.5+c*u+d. I am trying to solve it in Matlab but it does not find the solution. The problem is a body falling down. The body is under the effect of four forces: 1) weight (m*g) 2) flow in opposite direction of the weight (0.5*da*A*U0^2) 3) drag force (0.5*da*C_d*A*u^2), which would change verse depending on flow direction 4) acceleration (m*u') u(t=0)=0;
The C_d changes with velocity (relationship from the sphere in free flow C_d = 21.12/Reynolds+6.3/sqrt(Reynolds)+.25; % drag coefficient)
I coded this problem, but I can't get it to give a result, unless I fix C_d to a value. Would you be able to give a help, please?
da = 1.161; %density of air
g = 9.81; % Gravitational acceleration
d =1E-6*100; % characteristic length of the body
m = 1000*.957251*d.^3.09275; %Mass of the the body
A = 3.3108*d.^2.21672; %Surface area of the body
u_flow = .5; % Velocity air flow
syms u(t) c_D Re_ver
Re_ver = da*u(t)*d/1.8E-5; % Reynolds number
c_D = 21.12/Re_ver+6.3/sqrt(Re_ver)+.25; % drag coefficient;
a = m;
b = -.5*da*A*c_D;
c = m*g-.5*da*u_flow^2*A;
if c<0
c = -c;
end
u(t) = dsolve(a*diff(u,t) == b*u^2+c,u(0)==0);
t = linspace(0,10,500);
figure, plot(t,double(u(t)))
I am using Matlab R2014a
Thanks Antonio

8 个评论
显示 6更早的评论 隐藏 6更早的评论

John D'Errico
John D'Errico 2016-7-9
编辑：John D'Errico 2016-7-9
Why do you presume an analytical solution is there to be found in that case? You may need to use numerical tools.
Star Strider
Star Strider 2016-7-9
I would use the numeric ODE solvers, probably ode45, instead of the Symbolic math Toolbox.
mortain
mortain 2016-7-10
http://uk.mathworks.com/matlabcentral/answers/294607#comment_378361
Thanks John. I am not necessarily looking for an analytical solution. I tried int() but I could not get an answer. So I am interested to find a solution, even numerical.
http://uk.mathworks.com/matlabcentral/answers/294607#comment_378366 Hello Star, thanks for the reply. The problem with ode45 is that I am not capable to resolve it, maybe it's very simple and I am being stupid, might be a rookie question. I explicitly defined the coefficients of the equation but I get NaN:
if true
da = 1.161; %density of air
g = 9.81; % Gravitational acceleration
d =1E-6*100; % characteristic length of the body
m = 1000*.957251*d.^3.09275; %Mass of the the body
A = 3.3108*d.^2.21672; %Surface area of the body
u_flow = .5; % Velocity air flow
mu = 1.8E-5; %viscosity
%Re_ver = da*u*d/mu; % Reynolds number
%c_D = 21.12/Re_ver+6.3/sqrt(Re_ver)+.25; % drag coefficient;
%b = -.5*da*A*c_D;
a = m;
c = m*g-.5*da*u_flow^2*A;
if c<0
c = -c;
end
tspan = [0 5];
u0 = 0;
% Equation to solve a*u' = b*u^2 + c;
% I explicitly defined the coefficient as a constant or variable of u
[t,u] = ode45(@(t,u) ((-.5*da*A*(21.12/(da*u*d/mu)+6.3/sqrt(da*u*d/mu)+.25))*u^2+c)/a, tspan, u0);
plot(t,u,'-o')
end
mortain
mortain 2016-7-10
I solved it for now numerically. I found that, even though the equation is of the form: ...a/u*u^2... Matlab doesn't like the u at the denominator. By expanding the equation I could numerically solve it by doing a simple forward time marching. Please, let me have your thoughts, in case I am happy to further explain myself.
Walter Roberson
Walter Roberson 2016-7-11
I ran it through Maple, but I gave up on it after it had been computing the differential equation for several hours.
The Maple setup would be
Q := v->convert(v,rational);
da := Q(1.161);
g := Q(9.81);
d := 100*Q(0.1e-5);
m := 1000*Q(.957251)*d^Q(3.09275);
A := Q(3.3108)*d^Q(2.21672);
u_flow := Q(.5);
Re_ver := da*u(t)*d/Q(0.18e-4);
c_D := Q(21.12)/Re_ver+Q(6.3)/sqrt(Re_ver)+Q(.25);
a := m;
b := Q(-.5)*da*A*c_D;
c := m*g-Q(.5)*da*u_flow^2*A;
ut := dsolve({a*(diff(u(t), t)) = b*u(t)^2+c, u(0) = 0});
Raising values to fractional powers like 3.09275 often takes a lot of computation symbolically, as it would typically end up checking the 4000 complex roots of 3629/4000...
mortain
mortain 2016-7-12
Hello Walter, thanks a lot for your reply.
I do not have Maple, so I cannot test it. I am amazed that Matlab is not capable of solving this problem.
Assuming b = constant/u;
  1. If I'd formulate the problem in this way: u' = b*u^2; Matlab wouldn't be able to solve it
  2. If I'd formulate the problem in this way: u' = (constant)*u; Matlab would be able to solve it
I guess that it does not perform symbolic simplification before setting up a solution.
Thanks
Karan Gill
Karan Gill 2016-7-13
In the ODE a*u' = b*u^2+c , b has a term 1/sqrt(u). My guess is that this term complicates the calculation. The ODE is of the form:
C*u' = u^2(C + C/u + C/sqrt(C*u)) + C
|----'b' term-------|
where C represents the constants. I'm afraid I don't see how b = constant/u. Could you point out where I misunderstand?
Walter Roberson
Walter Roberson 2016-7-13
At your initial condition, u(0) = 0, your Reynold's number is 0. That is giving you a division by 0 for your drag coefficient, which is giving a singularity right from the start, which is preventing any resolution of the system.
I am not familiar with fluid dynamics, but when I look at various material it appears to me that the equation you are using might perhaps not apply at very low velocities.

请先登录，再进行评论。

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

回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

产品

标签

Community Treasure Hunt

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

Start Hunting!

Translated by