CODEWITHSUNDEEP
matlabode
user3648565
user3648565

Reputation: 79

Forcing matlab ODE solvers to use dy/dx = 0 IF dy/dx is negative

I need to numerically integrate the following system of ODEs:

dA/dR = f(R,A,B) dB/dR = g(R,A,B)

I'm solving the ODEs for a Initial-value stability problem. In this problem, the system is initially stable but goes unstable at some radius. However, whilst stable, I don't want the amplitude to decay away from the starting value (to O(10^-5) for example) as this is non-physical since the system's stability is limited to the background noise amplitude. The amplitude should remain at the starting value of 1 until the system destabilises. Hence, I want to overwrite the derivative estimate to zero whenever it is negative.

I've written some 4th order Runge-Kutta code that achieves this, but I'd much prefer to simply pass ODE45 (or any of the built in solvers) a parameter to make it overwrite the derivative whenever it is negative. Is this possible?

Upvotes: 2

Views: 396

Answers (1)

horchler
horchler

Reputation: 18484

A simple, fast, efficient way to implement this is via the max function. For example, if you want to make sure all of your derivatives remain non-negative, in your integration function:

function ydot = f(x,y)
ydot(1) = ...
ydot(2) = ...
...
ydot = max(ydot,0);

Note that this is not the same thing as the output states returned by ode45 remaining non-negative. The above should ensure that your state variables never decay.

Note, however, that that this effectively makes your integration function stiff. You might consider using a solver like ode15s instead, or at least confirming that the results are consistent with those from ode45. Alternatively, you could use a continuous sigmoid function, instead of the discontinuous, step-like max. This is partly a modeling decision.

Upvotes: 1

Related Questions