Estimate parameters of Ordinary Differential Equations (ODE)
23 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Jorge Armando Vazquez
am 19 Mär. 2021
Kommentiert: Jorge Armando Vazquez
am 19 Mär. 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 Kommentare
Akzeptierte Antwort
Alan Stevens
am 19 Mär. 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 Kommentare
Alan Stevens
am 19 Mär. 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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!