Solving ODEs in Matlab for a Non-System Problem

[Allamarein]

Let's suppose I know every coefficients this script:
and
y1 alfa
y2 v
y3 A
y4 T

function dy=isaacsimply(s,y)
dy = zeros(4,1)
global ...
dy(1)=(Cd_for*q_inf*h*(sin(y(1)))^2+...
g*y(3)*(rho-rho_inf)*cos(y(1))...
+E*U_inf*sin(y(3)))/...
(-rho*y(3)*y(2)^2); %dalfa
d_alfa=(Cd_for*q_inf*h*(sin(y(1)))^2+...
g*y(3)*(rho-rho_inf)*cos(y(1))...
+E*U_inf*sin(y(3)))/...
(-rho*y(3)*y(2)^2);

dy(2)=(q_inf*y(3)*sin(y(1))*cos(y(1))*d_alfa...
-g*y(3)*(rho-rho_inf)*sin(y(1))+...
E*U_inf*cos(y(1))-...
pi*h*rho*beta*(y(2)-U_inf*cos(y(1)))^2-...
y(2)*E)/(rho*y(2)*y(3)); %dv
d_v= (q_inf*y(3)*sin(y(1))*cos(y(1))*d_alfa...
-g*y(3)*(rho-rho_inf)*sin(y(1))+...
E*U_inf*cos(y(1))-...
pi*h*rho*beta*(y(2)-U_inf*cos(y(1)))^2-...
y(2)*E)/(rho*y(2)*y(3));

dy(3)=(E-y(3)*rho*d_v)/(rho*y(2)); %dA
d_A=(E-y(3)*rho*d_v)/(rho*y(2));

dy(4)=(E*Cp*y(4)+Nu/d*k*C*(T_inf-y(4))...
-rho*y(2)*Cp*y(4)*d_A-rho*y(3)*Cp*y(4)*d_v)...
/(rho*y(3)*Cp*y(2)); % dT

This script solves the system.jpg?

 

[Dr Transport]

I'd use Runge-Kutta on this one as long as you know the initial conditions...

 

[Allamarein]

I'd use Runge-Kutta on this one as long as you know the initial conditions...

ICs:

y1=pi/2;
y2=51.6;
y3=1.96e-7
y4=300

I'd know my transposition from mathematical problem to Matlab is correct. My results tell me I am wrong.

 

[Dr Transport]

At a quick look what you have listed as d_v is actually y(2)...

You also must look and see how well conditioned your functions are, i.e. your initial conditions are orders of magnitude different...That makes a huge difference in how well the solution is calculated.

 

[Allamarein]

At a quick look what you have listed as d_v is actually y(2)...

d_v=dy(2);
I used this expedient because I thought typing dy(2) directly, Matlab uses that dy(2) and not that calculated at precedent step. In a word it is a complicated story, but now I suppose defining d_v is useless.

You also must look and see how well conditioned your functions are, i.e. your initial conditions are orders of magnitude different...That makes a huge difference in how well the solution is calculated.

You'd be right, but these are my ICs. I supposed Matlab should give a warning (i.e an out tollerance value) if it finds troubles. If Matlab gives me a solution, I should trust that is correct (and fortunately I don't)
Anyway, maybe should I try to convert my system in a dimensionless one?

 

[Dr Transport]

Your equations are highly nonlinear, there are many issues with getting a robust solution. Linearize them to get a 1st order, then add in layers of complexity to ensure you don't get off track.

 

[Allamarein]

Your equations are highly nonlinear, [...]. Linearize them to get a 1st order, then add in layers of complexity to ensure you don't get off track.

It seem a good hints. May you suggest a first guess of a simpler system? How could I linearize my system.jpg? I am not very experienced in this kind of thing.
Anyway I am disapponited: I thought Matlab reports these kind of troubles when it solves ODE

 

[Dr Transport]

Since your initial conditions for $$\alpha$$ are at $$\pi/2$$ try expanding all of you trigonometric terms near that point to get an idea of how the set of equations behaves initially.


As for Matlab reporting troubles, it is just a numerical analysis tool, it doesn't know when a solution goes off in a wrong direction as long as it is getting numbers, it'll keep spitting them out until you tell it to stop. Judging the correctness of the solution is what you have to figure out. this is why I suggested working on linearizing first, you build up a confidence level in your solution space. You linearize, you set variables to the max and min that they will ever be to eliminate terms, etc... all of a sudden, you have mapped the solution space and you have understanding.

 

[Allamarein]

Since your initial conditions for $$\alpha$$ are at $$\pi/2$$ try expanding all of you trigonometric terms near that point to get an idea of how the set of equations behaves initially.

Dr Transport seems very stong in this kind of problems.
Please be patient with me: it' s the first time I tackle ODEs system and I should join the fray without a solid preparation.
Could you suggest a reference to learn how solve this kind of problem?

Anyway if I understand I should write a Taylor's expansion about $$\pi/2$$, therefore:
$$\sin(\alpha) \approx 1-1/2\, \left( \alpha-1/2\,\pi \right) ^{2}$$
$$\cos(\alpha) \approx -\alpha+1/2\,\pi$$
$$\sin ^2 (\alpha ) \approx 1 - \left( {\alpha - 1/2\,\pi } \right)^2$$
or better at the first order:
$$\sin (\alpha ) \approx 1 - {1 \over 2}\left( {{{\pi ^2 } \over 4} - \alpha \pi } \right)$$
$$\cos(\alpha) \approx -\alpha$$
and try to see where my plot goes.

If I understand my ODEs system is non-linear, so it has more than a solution. Hence Matlab could offer a solution different from that I can read on my reference (and that is confirmed in experimental tests).
I typed bunkum or more or less that's the point?
Let's suppose this is the point, nevertheless my ODEs is derived by governing equations. Therefore, given ICs, the behavior of my unknows (angle, velocity, area, temperature) should be only.
Probably I typed a lot of bunkum :-)

 

[f95toli]

No, non-linear does not mean that it has more than one solution.
It does however mean that you have to be very careful when you solve it because it will be very susceptible to numerical errors (and your system is likely to be chaotic).
You also need to be careful when you choose which solver the use, there is a reason for why there are so many different solvers available in Matlab, they are good a different things.

 

[Allamarein]

No, non-linear does not mean that it has more than one solution.
It does however mean that you have to be very careful when you solve it because it will be very susceptible to numerical errors (and your system is likely to be chaotic).
You also need to be careful when you choose which solver the use, there is a reason for why there are so many different solvers available in Matlab, they are good a different things.

I suspected this behavior. In fact I tried to use different solver (ode45, ode113, ode15s,...).
In order to learn ODE solving in Matlab I read a lot of tutorials, but they offer easy problems. I consider Matlab should be helpful when the system is complicated (like mine).
I wouldn't be a nuisance, but if it is only a numerical errors matter, I expect that Matlab gives me a warning when it has a problem, otherwise I never trust in its solutions.
Because in other situations, Matlab warns me when it has a large tollerance error and the solution si not reliable, I supposed it warns me in this case too, or they are two different circumstances?

Anyway I think to need a noble soul has a browse through my M.files or else a pithy tutorial.

 

[Dr Transport]

My suggestion is to expand the sinusiodal terms to low order, and go from there. Now, another trick to get some insight is to set all of the derivative terms on the RHS to a constant and plot the RHS as a function of $$\alpha$$ etc... You might see something interesting. Chaotic behavior in coupled equations comes about in very subtle ways, and given the complexity, I'd bet that numerical solution will take some time to nail down.

 

[f95toli]

In order to learn ODE solving in Matlab I read a lot of tutorials, but they offer easy problems. I consider Matlab should be helpful when the system is complicated (like mine).

But Matlab can never be more than a very good tool. It is a software package, not a mathematician-in-a-box.
The "easy problems" you refer to is a good start and teaches you how to use Matlab, definitely you need to understand how to solve "uncomplicated" problems before you tackle something as complicated as your system.
The main problem here is not the software; the problem is simply that the math/physics is very complicated.

 

[Allamarein]

I tried to use your wise hints.
I used taylor series for trigonometric functions
$$\eqalign{ & \sin \left( \alpha \right) \approx 1 - {1 \over 2}\left( {\alpha - {1 \over 2}\pi } \right)^2 \cr & \cos \left( \alpha \right) \approx - \alpha + {1 \over 2}\pi \cr}$$
I didn't discover a large difference.

I uploaded results for the case with taylor expansion and for the "exact" system.
The path.jpg is another story:
when I calculated $$\alpha$$ I can plot my path in X-Y plane, moving of a $$ds$$ in the direction suggested from the angle.
$$X\left( i \right) = X(i - 1) + \cos \left( \alpha \right)ds$$
$$Y\left( i \right) = Y\left( {i - 1} \right) + \sin \left( \alpha \right)ds$$

Black line is the path born by ODEs
Yellow,green and red paths are calculated with correlation equation or read from reference.
Black line should be similar at the red one ultimately.

 

[Dr Transport]

Now your getting there, keep plugging along, you've made progress.

 

[Allamarein]

Now your getting there, keep plugging along, you've made progress.

Plugging? Progress?
Sincerely I think to be in the middle of the ocean and I know where is the coast.
Would you like to see my M.files?

 

[Dr Transport]

Part of the learning/education process is working these issues to an end and becoming comfortable in working on your own. You think you have not gotten anywhere, but in all actuality, you have made a bunch of progress towards your goal.

 

[Allamarein]

Part of the learning/education process is working these issues to an end and becoming comfortable in working on your own. You think you have not gotten anywhere, but in all actuality, you have made a bunch of progress towards your goal.

You are right. I hope to achieve my goal quickly.
Early I submit you a new problem! It's very similar..but it is not a system.