Loop over ode45 to find minimum of a parameter

Question

Nader Mohamed 2021-10-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1574158-loop-over-ode45-to-find-minimum-of-a-parameter

评论： Star Strider 2021-10-28

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 个评论
显示无隐藏无

Walter Roberson 2021-10-28

What should happen if the entry with minimum errb is not the entry with the minimum errk ?

Nader Mohamed 2021-10-28

You're right. I edited the code to compute the error only compared to one analytical solution

My idea is then after having the vector err - find the min(errbk) and output the corrisponding b and k

请先登录，再进行评论。

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

Answer 1

Star Strider 2021-10-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1574158-loop-over-ode45-to-find-minimum-of-a-parameter#answer_819103

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 个评论
显示无隐藏无

Nader Mohamed 2021-10-28

Thanks! I understand that with lsqcurvefit I can fit the data with the results, but I still don't understand how can I retrieve b and k from that

Star Strider 2021-10-28

在 MATLAB Online 中打开

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.

.

请先登录，再进行评论。