Ordinary differential equation - pratical problem

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

[John D'Errico]

Why do you presume an analytical solution is there to be found in that case? You may need to use numerical tools.

[Star Strider]

I would use the numeric ODE solvers, probably ode45, instead of the Symbolic math Toolbox.

[mortain]

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]

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]

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]

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]

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]

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.