Estimate parameters of Ordinary Differential Equations (ODE)

23 次查看（过去 30 天）
Jorge Armando Vazquez
评论： 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);

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

采纳的回答

Alan Stevens
Alan Stevens 2021-3-19
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更早的评论
Alan Stevens
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by