Fitting with coupled differential equations

2017 12 20
1 个回答
回答已采纳
更新时间：2017 12 21
7 次查看（30 天）

0 个投票

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
Torsten 2017-12-20
Could you include a graphic of "pulses" over time ?
Best wishes
Torsten.
Silke
Silke 2017-12-20
Yes, sure. I have included a figure showing the full pulse and a zoom-in the interesting region.

请先登录，再进行评论。

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

 采纳的回答

Torsten
Torsten 2017-12-20

0 个投票

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
Silke 2017-12-20
Well, in principle pulses is a 1 * 3600 vector (as attached to the question at the top). I could make it an n*2 array though, if this helps. I will check the example "ODE with time-dependent terms". Thanks for the suggestion.
Silke
Silke 2017-12-20
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.
Torsten
Torsten 2017-12-20
编辑：Torsten 2017-12-20
If "pulses" were an (nx2)-matrix, you could use
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
Best wishes
Torsten.
Silke
Silke 2017-12-20
For some reason, I still get this error:
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.
Error while evaluating UIControl Callback
And the error occurs exactly the first time when it tries to evaluate
dCC(1) = pulse_actual*phi_n - k_1n .* x(1) - k_2 .* x(1) .* x(2);
I am not sure what this means. So the programm does not even try to evaluate dCC(2). Do you have any idea what might cause this error?
Torsten
Torsten 2017-12-20
How exactly does the function "DiffEqSolver" look like ?
Best wishes
Torsten.
Silke
Silke 2017-12-20
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
Torsten 2017-12-21
编辑：Torsten 2017-12-21
If you call
[t,x] = ode45(@(t,x) myode(t, x, pulses, phi_n, k_1n, k_2, k_1p), tspan, ic, opts);
the returned x will be a matrix of size (numel(tspan) x 2).
So calculating
f1x = 1.6e-19 * (mu_n * CCx(1)+mu_p*CCx(2));
makes no sense since CCx is a 2-dimensional matrix.
Best wishes
Torsten.
Silke
Silke 2017-12-21
I see your point. But what is confusing me is that the programm actually does not even get so far. The error occurs at the moment when it tries to evaluate this:
dCC(1) = pulse_actual*phi_n - k_1n .* x(1) - k_2 .* x(1) .* x(2);
I have added breakpoints to follow the programm and it does not even start evaluating dCC(2). If I remove the semicolon at the end, it also does not show any result in the command window. So I think there is an error earlier in the code.
Do you have any idea what might be wrong there?
Torsten
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
Silke 2017-12-21
Yes, indeed, this did the job. Thanks a lot for your help.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by