Re: lsode use

[Juan Pablo Carbajal]

Cöline,

Just a formal detail. Please answer at the bottom or between the lines
(as I am doing now). This makes it easy for readers to follow the
mailing list archives.


On Wed, Aug 26, 2015 at 10:36 AM, Céline <address@hidden> wrote:
> In fact, my J values are different at each time t and are calculated
> separately at each time.
> So I need, each time, to give my J value to my ODE to be able to resolve
> it...
>
> I've tried your advices, but I do probably something wrong since it does not
> work at all...
>
> Here is the code of the loop (I've all the initialization of parameters and
> variables before this):
>
> for t=0:step:ts
>      x=@(X0,J) humsolide(X0,t,J);
>     [X, ISTATE, MSG]=lsode(x,X0,temps);

What is temps?

>     psat=psat0*exp(-(Lm/Rg)*((1/T)-(1/T0)));
>     if X(i)>=Xc
>       %evaporation rate calculation
>       J=a*Mw*ke*((psat/(Rg*T))-(p/(Rg*T))*(Y/((Mw/Ma)+Y)));
>     elseif X(i)<=Xc
>       %parameters calculation
>       Ri=R*(((X-Xr)/(Xc-Xr))^(1/3));
>       ki=((eps*D)/(tau*((R)^2)))*(1/((1/Ri)-(1/R)));
>       %evaporation rate calculation
>       J=a*Mw*(1/((1/ke)+(1/ki)))*((psat/(Rg*T))-(p/(Rg*T))*(Y/((Mw/Ma)+Y)));
>   end
>   Y=((Ms*J)/G)+Y0;
>   T=T0-((Ms*J*Lk)/(cpa*G));
>   %Adding the results in "results matrix"
>   Yres=[Yres;Y];
>   Jres=[Jres;J];
>   i+=1;
> end
>
> And my function humsolide:
>
> function XDOT=humsolide(X0,t,J)
>   XDOT=-J;
>   endfunction
>
> Do you know what's wrong ?
>
>
Form what I understand you have a system of the form

dx/dt = -J(x)

With J switching functionla form baed on the values of x. Thisi is an
hybrid dynamical system and you can opt to solve it with lsode or with
the odepkg
I did not check how the to forms of J glue to each other (continuous?
1st derivative continuous?), but assuming J is smooth enough you can
use lsode in the following way

function J = Jvalue (x, params)
    # params is a struct with all your parameter values. Here comes a
block where you copy them ot variables of a nice name, if you want.
Otherwise use, e.g. parasm.Mw
    if x > params.Xc
        J = # function for x >= xc
    else # note that you elseif is not necessary and overlapping
        J = # function for x < xc
    endif
endfunction

function xdot = humsolide (x, t)
      xdot = Jvalue(x);
endfunction

Or you can integrate everything inside humsolide

function xdot = humsolide (x, t)
   # params is a struct with all your parameter values. Here comes a
block where you copy them ot variables of a nice name, if you want.
Otherwise use, e.g. parasm.Mw
    if x > params.Xc
        J = # function for x >= xc
    else # note that you elseif is not necessary and overlapping
        J = # function for x < xc
    endif
     xdot = J;
endfunction

and pass it to lsode as

func = @(x,t) humsolide (x,t, params)

The script will then be

#block where params is defined

func = @(x, t) humsolide(x,t,params);
[X, ISTATE, MSG] = lsode(x,X0, 0:step:ts);

As you see no need for the for loop.

If your J function is not smooth, then you should ask the question

Do I know the switching times a priori?

If the answer is yes, you can integrate using lsode and passing the
switching times as the argument T_CRIT of lsode (see documentation).

If the answer is no, then you should try to use odepkg, solver ode45,
and use the event function. Please refer to the documentation.

Good luck



>
> --
> View this message in context: 
> http://octave.1599824.n4.nabble.com/lsode-solver-problem-tp4672270p4672289.html
> Sent from the Octave - General mailing list archive at Nabble.com.
>
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://lists.gnu.org/mailman/listinfo/help-octave