Sensitivity analysis of an ODE

Question

Dereje 2019-1-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/439874-sensitivity-analysis-of-an-ode

评论： Dereje 2019-1-15

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

2 个评论
显示无隐藏无

madhan ravi 2019-1-15

Honestly you don't need a loop at all.

Dereje 2019-1-15

The loop was intended for a different purpose (Not related with my original question).

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jan 2019-1-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/439874-sensitivity-analysis-of-an-ode#answer_356502

编辑：Jan 2019-1-15

在 MATLAB Online 中打开

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).