Help with ODE system inside lsqnonlin

Hi All,

I am getting the following error when simulating a system of ODE inside lsqnonlin:

I have two ICs and it looks to me that the function I passed to ODE45 has two outputs and not 0 as the error says.

Below how I coded it.

Any help would be greatly appreciated!

Thank you in advance, Gabri

Fast_first3=fast(1:15121,:);
Slow_first3=y2(1:3,:)

options = optimoptions(@lsqnonlin, 'Algorithm', 'levenberg-marquardt','Display','final-detailed')
lb=[]; ub=[];
tic
[k]=lsqnonlin(@(k) estimation(k,Fast_first3,Slow_first3), k0, lb, ub, options);
toc

function Err=estimation(k, Fast_first3, Slow_first3)
global s0 s1 delta Ks R0 r1 gamma Rinf V0 vs dw 

%theta1_mine - r1 included
s0=k(1);
s1=k(2);
delta=k(3);
Ks=k(4);
R0=k(5);
r1=k(6);
gamma=k(7);
Rinf=k(8);
V0=k(9);

%Cycle A
tspanA_1 = linspace(0, (27.6-0.01), 27.6/0.01);
IC1_A=[10^-10 10^-10];
[t1_A, sol1_A] = ode45(@stateEqs1_A, tspanA_1, IC1_A); 

P=Ks.*(s0+s1.*sol1_A(:,1).^delta).*dw.*sol1_A(:,2);
Ra=(Rinf-(Rinf-R0).*exp(-sol1_A(:,1)./V0)).*(1+r1.*(pi.*dw.*sol1_A(:,2)./vs).^gamma);
Dw=(2.*sol1_A(:,1))./(pi*dw);

%P(5061)=[]; P(10121)=[]; P(15181)=[]; %Remove last measurement in fast
Ra=[Ra(5061,:); Ra(10122,:); Ra(15183,:)];
Dw=[Dw(5061,:); Dw(10122,:); Dw(15183,:)];

 E1=sqrt(inv(10^4)).*(Fast_first3-P);
    E2=sqrt(inv(10^-4)).*(Slow_first3(:,1)-Ra);
    E3=sqrt(inv(10^-4)).*(Slow_first3(:,3)-Dw);

 Err=[E1; E2; E3];
 
 function output = stateEqs1_A(t, X)
%global G1 g s0 s1 delta Ks R0 r1 gamma Rinf V0 vs dw ds0
u=.0106;

x1	=	X(1)	;
x2	=	X(2)	;

output=[pi*dw*x2;
    (vs*(u-x2))/(dw*(s0+s1*x1^(delta)))];


 end
 
end

[Jan]

sqrt(inv(10^-4)) ?! What about writing 100 instead of wasting time for INV and SQRT?

Remember that 10^-10 is an expensive power operation, while 1e-10 is a cheap constant.

[Gabriele Galli]

Thank you for the advice! I copied them from another file I coded using matrices and forgot to change that twisted sintax. Also, your advice in the answer was really helpful. The issue was related to dw. Thank you very much again! Gabri

Answer 1

[Jan]

It depends on what dw and the other global variables are.

Set a break point in this line:

output = [pi * dw * x2; ...
         vs * (u - x2) / (dw * (s0 + s1 * x1^delta))];

and check the values of the variables.

Providing dw as global variable and then re-using it in a nested function is prone to bugs.