In the following example, we solve the Ordinary Differential Equation $$
y''-\mu(1-y^2)y'+y=0,
$$
with the initial condition $y(0)=$, $y'(0)=0$, and $\mu=1$. Since the original equation is of second order we define the state of the equivalent first order equation as $v=(y,y')$. The right hand side is computed by the function %SUN_vdp1
(defined in the Sundials module):
function vdot=%SUN_vdp1(t, v)
mu = 1;
vdot = [v(2); mu*(1-v(1)^2)*v(2)-v(1)]
end
[t,v] = cvode(%SUN_vdp1, [0 20], [0; 2]);
plot(t',v')
CVODE can solve initial value problems of the kind
$$
\begin{array}{rcl}y^\prime & = & f(t,y,p), \quad y(t_0) =y_0(p),\\
s^\prime & = & \frac{\partial f}{\partial y}s+\frac{\partial f}{\partial p},
\quad s(t_0) = \frac{\partial y_0(p)}{\partial p},\end{array}
$$
where the sensitivity function is $s={\partial y}/{\partial p}$. The right hand side of the sensitivity equation can be provided
by the user (see the sensRhs
option in the Sensitiviy
help page) but by default the directional derivatives are approximated by finite differences. The call to cvode
takes the form
cvode(fun, tspan, y0, sensPar = p)
When using a Scilab function, fun
has the prototype dydt = fun(t,y,p,...)
. Here is an example where we compute the sensitivity of the solution of the Van Der Pol equation with respect to parameter $\mu$:
function vdot=vdp(t, y, μ)
vdot = [y(2); μ*(1-y(1)^2)*y(2)-y(1)]
end
[t,y,s] = cvode(vdp, 0:0.1:10, [2;1], sensPar = 1);
plot(t',s')
legend("$\Large\partial y_"+["1","2"]+"/\partial μ$");