Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
Drexel University or The Math Forum.
|
|
Torsten
Posts:
1,182
Registered:
11/8/10
|
|
Re: Stiffness problem in ODEs
Posted:
Mar 6, 2013 2:43 AM
|
|
"Fotis Dimitrakopoulos" <f.dimitrak@gmail.com> wrote in message <kh4ltv$8ub$1@newscl01ah.mathworks.com>...
> Hello every one, I am trying to numerically solve a system of 5 coupled diff. equations using ode15s (i have also tried using different functions), but I get the following warning:
>
> "Warning: Failure at t=2.289553e-01. Unable to meet integration tolerances without reducing the step
> size below the smallest value allowed (4.440892e-16) at time t. ?
>
> Could anyone have an idea about what could I do to get a result? My code is the following:
>
> %define functions
> classdef beta
>
> methods(Static)
> function y = g1(gp, gs, g, yt)
> y = 41*gp^3/(96*pi^2) + (1/16*pi^2)^2*((44/3)*gs^2*gp^3 + (9/2*g^2 + (199/18)*gp^2 - (17/6)*yt^2)*gp^3);
> end
> function y = g2(gp, gs, g, yt)
> y = -19*g^3/(96*pi^2) + (1/16*pi^2)^2*(12*gs^2*g^3+((35/6)*g^2 + (3/2)*gp^2 - (3/2)*yt^2)*g^3);
> end
> function y = g3(gp, gs, g, yt)
> y = -7*gs^3/(16*pi^2) + (1/16*pi^2)^2*gs^3*((9/2)*g^2 - 26*gs^2 + (11/6)*gp^2 - 2*yt^2);
> end
> function y = gamma(g, gp, gs, yt, lam)
> y = (-9*g^2 - 3*gp^2 + 12*yt^2)/(64*pi^2) + ((1/(16*pi^2))^2)*((1/6)*lam^2 ...
> - (27/4)*yt^4 + 20*gs^2*yt^2 + (45/8)*g^2*yt^2 + (85/24)*gp^2*yt^2 - (271/32)*g^4 + (9/16)*g^2*gp^2 + (431/96)*gp^4);
> end
> function y = bf(yt, g, gp)
> y = 3/16*(3*g^4 + 2*g^2*gp^2 + gp^4) - 3*yt^4;
> end
> function y = blam(yt, g, gp, gs, lam)
> y = 4*lam*beta.gamma(g, gp, gs, yt, lam) + (12*lam^2 + beta.bf(yt, g, gp))/(8*pi^2) + (1/(16*pi^2))^2*(80*gs^2*yt^2*lam - 192*gs^2*yt^4 + (915/8)*g^6 - (289/8)*gp^2*g^4 - (27/2)*yt^2*g^4 - (73/8)*lam*g^4 - (559/8)*gp^4*g^2 + 63*gp^2*yt^2*g^2 + (39/4)*gp^2*lam*g^2 - 3*yt^4*lam + (45/2)*yt^2*lam*g^2 - (379/8)*gp^6 + 180*yt^6 - 16*gp^2*yt^4 - (26/3)*lam^3 - (57/2)*gp^4*yt^2 - 24*yt^2*lam^2 + 6*(3*g^2 + gp^2)*lam^2 + (629/24)*gp^4*lam + (85/6)*gp^2*yt^2*lam);
> end
> function y = yuk(yt, g, gp, gs, lam)
> y = yt*(9/2* yt^2 - 9*g^2/4 - 17/12*gp^2 - 8*gs^2)/(16*pi^2) + ((1/16*pi^2)^2)*yt*(-108*gs^4 + 9*g^2*gs^2 + (19/9)*gp^2*gs^2 + 36*yt^2*gs^2 - (3/4)*gp^2*g^2 - (23/4)*g^4 + (1187/216)*gp^4 - ...
> 12*yt^4 + (1/6)*lam^2 + yt^2*((225/16)*g^2 + (131/16)*gp^2 - 2*lam));
> end
> end
> end
>
>
> %system to solve
> function dy = system(t,y)
>
> dy = zeros(5,1);
> dy(1)=beta.g1(y(1),y(3),y(2),y(4));
> dy(2)=beta.g2(y(1),y(3),y(2),y(4));
> dy(3)=beta.g3(y(1),y(3),y(2),y(4));
> dy(4)=beta.yuk(y(4),y(2),y(1),y(3),y(5));
> dy(5)=beta.blam(y(4),y(2),y(1),y(3),y(5));
> end
>
> and then I give the command:
>
> [T,Y] = ode15s(@system,[0. 12],[0.357452 0.65171 1.1645 0.967834 0.139429]);
>
> P.S: This is my first post/question, so I apologize for any inconvenience.
>
> Thanks a lot
Your system of ODEs runs without problems with a different numerical integrator
on my computer.
Maybe the structure out of which you get the derivatives "dy" causes problems.
I'd suggest you test whether directly prescribing the derivatives without calling the beta.xx - functions works.
Best wishes
Torsten.
|
|
|
|
