MATLAB: Sensitivity analysis of an ODE

sensitivity analysis

Hi,
I have been asked to do the sensitivity analysis of an ode Model. The idea is new to me and I need some help to give me some hints.
The code is as follows: 
zspan=[0,400];
v0mat = [1 0.01 1];
zsol = {};
v1sol = {};
v2sol = {};
v3sol = {};
for k=1:size(v0mat,1)
v0=v0mat(k,:);
[z,v]=ode45(@rhs,zspan,v0);
zsol{k}=z;
v1sol{k}=v(:,1);
v2sol{k}=v(:,2);
v3sol{k}=v(:,3);
end
for r=1:length(v2sol)
q(r)=r;
end
for k1 = 1:length(v2sol)
  zsol04(k1) = interp1(v2sol{k1}, zsol{k1}, 0.4);
end
function parameters=rhs(z,v)
          alpha=0.08116; 
          db= 2*alpha-(v(1).*v(3))./(2*v(2).^2);
          dw= (v(3)./v(2))-(2*alpha*v(2)./v(1));
          dgmark= -(2*alpha*v(3)./v(1)); 
 parameters=[db;dw;dgmark];
end

Best Answer

An analysis of the sensitivity can be achieve by two methods:
1. Vary the initial values y0 and/or parameters p by a "small" amount delta. Calculate the difference of the results. The senisitivity is the quotient:
(y(end) - y_varied(end)) ./ delta
The size is criticial, because if it is too large, you get a too rought discretization, and if it is too small, you get mainly the cancellation error, which is amplified by the division by a tiny number. The goal is to get a y-y_varied in the magnitude of sqrt(eps). This is called "external numerical differentiation", because you measure the gradient externally.
2. Modify the integrator internally: In each step of the solver create the matrix by the numerical differentiation. Vary the parameters and current positions. You get a matrix and if you start with the unit matrix and multiply these sensitivity matrices cumulatively, you get a more accurate sensitivity matrix.
The methods can answer substantially different results, if there are poles or discontinuities near to the trajectory. For smooth functions the resulting matrices are identical in first order.
By the way, if you calculate the "Wronski" matrix by the internal numerical differentiation, this valuable information could be used in the step size controller also. The implemented controller in e.g. ODE45 varies the order of the integration scheme to measure the deviation of the result in each step. This is very similar to varying the step size or the positions or parameters to determine the sensitivity.
I assume, this is a homework and the 1st method is sufficient already. You have to run the integration with v0 and then 3 times with a varied v0(k) with k=1:3. Maybe you want to vary alpha also. Then use Answers: Anonymous functions for parameters (link).
Related Question