Estimate parameters of Ordinary Differential Equations (ODE)

Question

Jorge Armando Vazquez 2021-3-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/777962-estimate-parameters-of-ordinary-differential-equations-ode

评论： Jorge Armando Vazquez 2021-3-19

Hello, this code was from a post of some years ago. I am trying to estimate parameters of Ordinary Differential Equations (ODE). But I am receiveing the following error:

Subscripted assignment dimension mismatch.

Error in fminsearch (line 190)

fv(:,1) = funfcn(x,varargin{:});

I can´t solve it and I really need help

..

% == ODE ==

function dx=LV(t,x,theta)

dx=zeros(2,1);

alpha=theta(1);

beta=theta(2);

gamma=theta(3);

delta=theta(4);

dx(1) = x(1)*(alpha-beta*x(2));

dx(2) = -x(2)*(gamma-delta*x(1));

end

% == Error function ==

function err = ODE_fit(exp_t, exp_y, theta)

% exp_y = Experimental observation at time exp_t

[t,y] = ode45(@(t,X)LV(t,X,theta), exp_t, [5 3]);

err = sum((y-exp_y).^2); % compute error between experimental y and fitted y

end

% == Script ==

exp_t = [0, 0.20, 0.400, 0.60, 0.80, 1, 1.20, 1.40, 1.60, 1.80, 2];

exp_y = [4.35,4.26,2.96,3.13,2.25,2.65,3.22,2.85,4.97,4.99,5.94;...

2.60,3.23,2.50,1.89,2.25,1.21,1.05,1.55,1.66,1.44,1.76]';

theta0 = [2.5 0.75 4.25 1.5];

p_estimate = fminsearch(@(theta)ODE_fit(exp_t, exp_y, theta), theta0);

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Alan Stevens 2021-3-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/777962-estimate-parameters-of-ordinary-differential-equations-ode#answer_652942

在 MATLAB Online 中打开

This works (see the comments):

% == Script ==
exp_t = [0, 0.20, 0.400, 0.60, 0.80, 1, 1.20, 1.40, 1.60, 1.80, 2];
exp_y = [4.35,4.26,2.96,3.13,2.25,2.65,3.22,2.85,4.97,4.99,5.94;...
         2.60,3.23,2.50,1.89,2.25,1.21,1.05,1.55,1.66,1.44,1.76]';
theta0 = [2.5 0.75 4.25 1.5];
p_estimate = fminsearch(@(theta)ODE_fit(theta, exp_t, exp_y), theta0);
disp(p_estimate)
% == Error function ==
function err = ODE_fit(theta, exp_t, exp_y)  %%%%%% Must have theta first
% exp_y = Experimental observation at time exp_t
[~,y] = ode45(@(t,X)LV(t,X,theta), exp_t, [5 3]); 
err = sum((y-exp_y).^2,'all'); %%%%%% Use 'all' to sum over rows and columns
end
% == ODE ==
function dx=LV(~,x,theta)
dx=zeros(2,1);
alpha=theta(1); 
beta=theta(2); 
gamma=theta(3); 
delta=theta(4);
dx(1) = x(1)*(alpha-beta*x(2));
dx(2) = -x(2)*(gamma-delta*x(1));
end

3 个评论
显示 1更早的评论隐藏 1更早的评论

Alan Stevens 2021-3-19

Not for me it doesn't! How are you running it? Copy and paste into a new script file, save it and then run it. The output I get from disp(p_estimate) is

2.4346 1.1923 2.5901 0.6243

Jorge Armando Vazquez 2021-3-19