ode45 + lsqcurvefit- multiple ODEs
3 views (last 30 days)
Show older comments
Camila de Paula
on 18 Jun 2018
Answered: Star Strider
on 18 Jun 2018
Hi,
I am trying to fit my experimental data to a model I have created that include 2 coupled ODEs. I created a separate function for defining the 2 ODEs, where a is the array that contains the 2 constants I'm trying to fit:
function dy = CalcCoupODEs(t,y,a) %dy= matrix with all 2 ODEs
dy = zeros(2,1); % variable definition
dy(1) = 1.41.*a(1) .* ((y(2)*y(1)^2)^2)*y(1) - a(2) .* y(1)*t ; %ODE for R
dy(2) = -a(1)*(y(2)*y(1)^2)^2; %ODE for n
end
Then, I call the function in the following way:
cycle_DMACl = [20, 25, 30, 40, 45]; % this is t variable
R_DMACl = [18.34, 11.29, 7.09, 6.51, 4.396]; %this is y(1)
n_DMACl = [92.31, 61.9, 53.82, 26.04, 11.87]; % this is y(2)
a0 = [1, 1]; %initial guess for constants
lb1 = [0, 0]; %lower bound
a = lsqcurvefit(@CalcCoupODEs, a0, cycle_DMACl(:), R_DMACl(:), n_DMACl(:), lb1);
I get the following error: Warning: Length of lower bounds is > length(x); ignoring extra bounds.
It's obviously not trying to fit the correct thing (a and lb are both 1x2)
What am I doing wrong and what is the b est way to do this?
Thank you!
0 Comments
Accepted Answer
Star Strider
on 18 Jun 2018
cycle_DMACl = [20, 25, 30, 40, 45]; % this is t variable
R_DMACl = [18.34, 11.29, 7.09, 6.51, 4.396]; %this is y(1)
n_DMACl = [92.31, 61.9, 53.82, 26.04, 11.87]; % this is y(2)
Rn_DMACl = [R_DMACl(:), n_DMACl(:)]; % Create Single Matrix
a0 = [1, 1]; %initial guess for constants
lb1 = [0, 0]; %lower bound
a = lsqcurvefit(@(b,t)objfcn(b,t), [a0 0 0], cycle_DMACl(:), Rn_DMACl)
function yest = objfcn(b,t)
a(1) = b(1);
a(2) = b(2);
ic = [b(3); b(4)];
function dy = CalcCoupODEs(t,y) %dy= matrix with all 2 ODEs
dy = zeros(2,1); % variable definition
dy(1) = 1.41.*a(1) .* ((y(2).*y(1).^2).^2).*y(1) - a(2) .* y(1).*t; %ODE for R
dy(2) = -a(1).*(y(2).*y(1).^2).^2; %ODE for n
end
[t,yest] = ode15s(@CalcCoupODEs, t, ic);
end
Rn_fit = objfcn(a,cycle_DMACl(:))
figure(1)
plot(cycle_DMACl(:), Rn_DMACl, 'p')
hold on
plot(cycle_DMACl(:), Rn_fit, '-r')
hold off
grid
The problem is that your ‘CalcCoupODEs’ keeps crashing, with the error:
Warning: Failure at t=2.000000e+01. Unable to meet integration tolerances without reducing the step size below
the smallest value allowed (5.684342e-14) at time t.
So it is encountering a singularity, and the integration stops. You need to fix that. I do not know what system of differential equations you are coding (you did not post them).
My code should work otherwise, although I cannot be certain of that, since the ODE integration crashes on the first time step, even with non-zero and non-identical initial conditions.
0 Comments
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!