Loop over ode45 to find minimum of a parameter
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
I'm trying to loop over an ode45 for different b and k, to find the couple of the two that minimize the error from the analytical solution. But when I run this code it enters in an infinite loop. What am I doing wrong?
T = readtable('samples.csv'); % three column [time,analytical_sol1,analytical_sol2]
test = @(t,y,b,k) [0;0;k/J1 * y(2) - k/J1 * y(1); T0/J2 - b*y(4)/J2 - (k/J2)*(y(2)-y(1))]; %my ode
err = []; % initialize error vector
for b = 0:0.01:10 %loop over different b
for k = 0:0.1:100 %loop over different k
[t,y] = ode45(@(t,y) test(t,y,b,k) , [0:0.01:10] , [0,0,0,0]); %solve my ode
errbk = abs( norm( T{:,3} - y(:,4) ) ); % compute error from the analytical solution 1
err = [err;b,k,errbk];
end
end
% then i would find b and k with the minimum errb and errk
2 Kommentare
Walter Roberson
am 28 Okt. 2021
What should happen if the entry with minimum errb is not the entry with the minimum errk ?
Antworten (1)
Star Strider
am 28 Okt. 2021
There are several examples on fitting differential equations to data, one being Coefficient estimation for a system of coupled ODEs — not trivial, however also not difficult.
.
2 Kommentare
Star Strider
am 28 Okt. 2021
It would be relatively straightforward to adapt my code to calculate ‘b’ and ‘k’. They become parameters, so the ‘kinetics’ function becomes —
function C=kinetics(theta,t,T0,J2)
% c0=[1;0;0;0];
c0 = theta(3:6);
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c) % k = theta(1), b = theta(2)
dcdt=zeros(4,1);
dcdt(1)= 0;
dcdt(2)= 0;
dcdt(3)= theta(1)/J1 * y(2) - theta(1)/J1 * y(1);
dcdt(4)= T0/J2 - theta(2)*y(4)/J2 - (theta(1)/J2)*(y(2)-y(1));
dC=dcdt;
end
C=Cv;
end
This is my best guess on how to implement your system of differential equations with my existing code. Here, the parameter vector ‘theta’ has ‘k’ and ‘b’ as the first two elements, and the initial conditions for the system of differential equations as the last four elements. All will be estimated by the optimisation funciton (lsqcurvefit, ga, or others). It may be necessary to edit this, because I do not understand what the objective is.
The ‘C’ output will be the result that matches the data to be regressed against. All will be matrices of column vectors.
.
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!