Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Solving ODEs system at Matlab

  1. Nov 11, 2008 #1
    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?
     

    Attached Files:

  2. jcsd
  3. Nov 11, 2008 #2

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    I'd use Runge-Kutta on this one as long as you know the initial conditions......
     
  4. Nov 12, 2008 #3
    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.
     
  5. Nov 13, 2008 #4

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  6. Nov 14, 2008 #5
    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'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?
     
    Last edited: Nov 14, 2008
  7. Nov 14, 2008 #6

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  8. Nov 15, 2008 #7
    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
     
  9. Nov 15, 2008 #8

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    Since your initial conditions for [tex] \alpha [/tex] are at [tex] \pi/2[/tex] 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.
     
  10. Nov 15, 2008 #9
    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 [tex] \pi/2 [/tex], therefore:
    [tex] \sin(\alpha) \approx 1-1/2\, \left( \alpha-1/2\,\pi \right) ^{2} [/tex]
    [tex] \cos(\alpha) \approx -\alpha+1/2\,\pi [/tex]
    [tex] \sin ^2 (\alpha ) \approx 1 - \left( {\alpha - 1/2\,\pi } \right)^2 [/tex]
    or better at the first order:
    [tex]\sin (\alpha ) \approx 1 - {1 \over 2}\left( {{{\pi ^2 } \over 4} - \alpha \pi } \right) [/tex]
    [tex] \cos(\alpha) \approx -\alpha [/tex]
    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 :-)
     
    Last edited: Nov 15, 2008
  11. Nov 15, 2008 #10

    f95toli

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  12. Nov 16, 2008 #11
    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.
     
    Last edited: Nov 16, 2008
  13. Nov 16, 2008 #12

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    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 [tex] \alpha [/tex] 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.
     
  14. Nov 16, 2008 #13

    f95toli

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  15. Nov 16, 2008 #14
    I tried to use your wise hints.
    I used taylor series for trigonometric functions
    [tex]
    \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}
    [/tex]
    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 [tex] \alpha [/tex] I can plot my path in X-Y plane, moving of a [tex] ds [/tex] in the direction suggested from the angle.
    [tex] X\left( i \right) = X(i - 1) + \cos \left( \alpha \right)ds [/tex]
    [tex] Y\left( i \right) = Y\left( {i - 1} \right) + \sin \left( \alpha \right)ds [/tex]

    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.
     

    Attached Files:

  16. Nov 16, 2008 #15

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    Now your getting there, keep plugging along, you've made progress.
     
  17. Nov 16, 2008 #16
    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?
     
  18. Nov 16, 2008 #17

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  19. Nov 17, 2008 #18
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Solving ODEs system at Matlab
Loading...