Estimate parameters of Ordinary Differential Equations (ODE)

조회 수: 46 (최근 30일)
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
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
Jorge Armando Vazquez
Jorge Armando Vazquez 2021년 3월 19일
You are right!!! Thank you very much!!

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Quadratic Programming and Cone Programming에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by