Estimate parameters of Ordinary Differential Equations (ODE)
23 views (last 30 days)
Show older comments
Jorge Armando Vazquez
on 19 Mar 2021
Commented: Jorge Armando Vazquez
on 19 Mar 2021
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 Comments
Accepted Answer
Alan Stevens
on 19 Mar 2021
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 Comments
Alan Stevens
on 19 Mar 2021
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
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!