Fitting with coupled differential equations

Question

Silke 2017-12-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/373873-fitting-with-coupled-differential-equations

评论： Silke 2017-12-21

采纳的回答： Torsten

pulses.mat

在 MATLAB Online 中打开

Hi there,

I am trying to fit experimental data with a set of coupled differential equations (CDE). At the moment, using these CDEs does not give an error message, but also does not fit my data.

function dCC = DiffEq(t, x);
    xdot(1,:) = phi_n - k_1n * x(1) - k_2 * x(1) * x(2);
    xdot(2,:) = phi_n - k_1p * x(2) - k_2 * x(2) * x(1);
    dCC = xdot;
end
[T,CCx] = ode45(@DiffEq, [t(1): dx : t(end)], [0 0]);
[xfit,resnorm, Jacob, CovB, MSE] = nlinfit(  handles.timecorr,handles.datacorr',@DiffEqSolver, handles.x0 );

The problem I am having now, is that in principle the first term of both CDEs should be only non-zero in a certain time window, which is defined by pulses. If I just multiply phi_n with pulses in the first term of the differential equation, I end up with a matrix, which is not allowed in ode45. Is there any way to overcome this issue?

Thank you!

2 个评论
显示无隐藏无

Torsten 2017-12-20

Could you include a graphic of "pulses" over time ?

Best wishes

Torsten.

Silke 2017-12-20

Yes, sure. I have included a figure showing the full pulse and a zoom-in the interesting region.

请先登录，再进行评论。

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

Answer 1

Torsten 2017-12-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/373873-fitting-with-coupled-differential-equations#answer_297109

So I suspect that "pulses" is an (nx2) array with the first column being "time" and the second column being values between 0 and 1 ?

If this is the case, take a look at the example

"ODE with time-dependent terms"

under

https://de.mathworks.com/help/matlab/ref/ode45.html

It can directly be applied to your ODE-system.

Best wishes

Torsten.

10 个评论
显示 8更早的评论隐藏 8更早的评论

Silke 2017-12-20

在 MATLAB Online 中打开

I am not sure what I am doing wrong, but I get now an error message saying:

Error using nlinfit (line 205) Error evaluating model function 'DiffEqSolver'.

Error in TRMC_GuiFit_SD_DiffEqfit_1>FITANDPLOT_Callback (line 404) [xfit,resnorm, Jacob, CovB, MSE] = nlinfit( handles.timecorr,handles.datacorr',@DiffEqSolver, handles.x0 ); %, handles.lb, handles.ub

Error in gui_mainfcn (line 95) feval(varargin{:});

Error in TRMC_GuiFit_SD_DiffEqfit_1 (line 42) gui_mainfcn(gui_State, varargin{:});

Error in matlab.graphics.internal.figfile.FigFile/read>@(hObject,eventdata)TRMC_GuiFit_SD_DiffEqfit_1('FITANDPLOT_Callback',hObject,eventdata,guidata(hObject))

Caused by: Index exceeds matrix dimensions.

And this is the code:

function dCC = myode(t, x, pulses, phi_n, k_1n, k_2, k_1p)
        xdot(1,:) = pulses.*phi_n - k_1n .* x(1) - k_2 .* x(1) .* x(2);
        xdot(2,:) = pulses.*phi_n - k_1p .* x(2) - k_2 .* x(2) .* x(1);
 end
tspan = [t(1): dx : t(end)];
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
ic = 0;
[t,x] = ode45(@(t,x) myode(t, x, pulses, phi_n, k_1n, k_2, k_1p), tspan, ic, opts);

pulses is still the same vector as included in the first question.

Silke 2017-12-20

在 MATLAB Online 中打开

function dCCconv = DiffEqSolver(param, t);
RC = 18e-9;
mu_n = param(5);
mu_p = param(6);
k_1n = param(1);
k_2 = param(2);
k_1p = param(3);
phi_n = param(4);
   filename1 = 'pulse300.dat';
   pathname1 = 'D:\My Documents\Experiments\TRMC\PulseShapes\';
   pulse = importfile1([pathname1, filename1]);
[P,I] = max(pulse(:,1)); %find the coordinates of the maximum value of t
[P1, I1] = find(pulse(:,2)>=0);
pulse_begin_time = pulse(I1(1),1);
pulse_begin = min(find(t>=pulse_begin_time));
pulse_end = max(find(t<=P));
[P3, I3] = max(t);
dx = t(2)-t(1);
t_x = t(pulse_begin : pulse_end);
pulses(:,2)=spline(pulse(:,1), pulse(:,2), t_x); % make the pulse and the experimental data have the same time scale
pulses(end+1:I3,2)= 0;
pulses(:,1) = t;
function dCC = myode(t, x, pulses, phi_n, k_1n, k_2, k_1p ) 
    dCC = zeros(2,1 ); 
pulse_actual = interp1(pulses(:,1),pulses(:,2),t ); 
dCC(1) = pulse_actual*phi_n - k_1n .* x(1) - k_2 .* x(1) .* x(2); 
dCC(2) = pulse_actual*phi_n - k_1p .* x(2) - k_2 .* x(2) .* x(1); 
end
tspan = [t(1): dx : t(end)];
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
ic = 0;
[t,x] = ode45(@(t,x) myode(t, x, pulses, phi_n, k_1n, k_2, k_1p), tspan, ic, opts);
CCx = x;
f1x = 1.6e-19 * (mu_n * CCx(1)+mu_p*CCx(2));
f2x = conv(exp(-t/RC),f1x); %convolution to take into acount of the time response of the cavity cell
f(1:I3) = f2x(1:I3)/72.501; %72.501 is taken from igor, maximum value of the integration of the cavity response
dCCconv = f;
assignin('base', 'f',f);
end

Torsten 2017-12-21

编辑：Torsten 2017-12-21

"ic" must be a vector with two components since you solve two differential equations. Thus "ic=0" must be replaced by something like ic=[0;0].

Best wishes

Torsten.

Silke 2017-12-21

Yes, indeed, this did the job. Thanks a lot for your help.

请先登录，再进行评论。

Fitting with coupled differential equations

2 个评论
显示无隐藏无

采纳的回答

10 个评论
显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Fitting with coupled differential equations

2 个评论 显示 无隐藏 无

采纳的回答

10 个评论 显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

10 个评论
显示 8更早的评论隐藏 8更早的评论