Re: Want to set dependent variable in system of ODEs being solved by ode45 to zero

[c.]

On 4 Nov 2012, at 06:57, clustro wrote:

> Hi Octave forum,
> 
> I have a problem with a system of differential equations I am trying to
> solve.
> 
> It is a system of nonlinear equations, discretized over length, and solved
> forwards in time.
> 
> The problem is, when another variable crosses a certain value, I want the
> dependent variable of a certain ODE to become 0 immediately and forever.
> 
> Unfortunately, I can't get ode45 in Octave to do it; it appears that
> internally in the ODE45 code, the previous values for the integrated
> variables are being remembered, so putting the zero-condition in my ODE
> doesn't do anything at all.
> 
> Observe code:
> 
> function [t,Cp] = pmin
> 
> r = [0:0.5:10]';
> Cp0 = 10*exp(-0.1*r);
> kI = 0.1;
> deltap = 0.1;
> kr = 0.1;
> params.r = r;
> params.kI = 1;
> params.deltap = 0.1;
> params.kr = 0.1;
> [t,Cp] = ode45(@(t,Cp) pCpODE(t,Cp,params),[0:5:250]',Cp0);
> [tt,rr] = meshgrid(t,r);
> mesh(tt,rr,Cp')
> 
> function dCpdt = pCpODE(t,Cp,params)
> r = params.r;
> kI = params.kI;
> deltap = params.deltap;
> kr = params.kr;
> n = length(Cp);
> dCpdt = zeros(n,1);
> R = r(end) - kr*t; <- R changes with time; defines zero condition below
> for i = 1:n
>       if r(i) > R
>               Cp(i) = 0; <- Right here is where I set the dependent variable 
> to zero
>       end
>       dCpdt(i) = -kI*Cp(i).*exp(-(R - r(i))/deltap);
> end
> 
> But it doesn't work. I don't know how to post graphs here, but clearly the
> result is not being zeroed. I thought that making Cp global or persistent
> would help, but it didn't do anything. Does anyone know how to get that to
> work? Or perhaps coding my own RK4 routine is the way to do it, e.g. so I
> can make sure the dependent variable gets set to zero in the solver itself?
> 
> Thanks!!


There seem to be misunderstandings at two different levels here

1) Mathematically tour problem is NOT an Ordinary Differential Equation, so 
ODE45 is most likely not the best tool to solve it. Maybe you can reformulate 
it as a Differential Algebraic Equation and solve it with DASPK, but I'm not 
sure ...

2) In the matlab/Octave interpreted language ALL functions have pass-by-vaule 
semantics, so changing the value of an input parameter inside a funcion
   leaves no effect after the function has returned

function b = fun (a)
  a = a + 1;
  b = a;
  printf ("in fun, a = %g", a)
endfunction

a = 1;
b = fun (a);

printf ("after running fun, a = %g", a)

HTH,
c.