Fitting with coupled differential equations

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!

[Torsten]

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

Best wishes

Torsten.

[Silke]

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

Answer 1

[Torsten]

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.

[Silke]

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]

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]

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]

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]

How exactly does the function "DiffEqSolver" look like ?

Best wishes

Torsten.

[Silke]

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]

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]

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]

"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]

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