Help with ODE system inside lsqnonlin

3 次查看（过去 30 天）
Gabriele Galli
Gabriele Galli 2021-6-25
评论： Gabriele Galli 2021-6-25
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
  2 个评论
Jan
Jan 2021-6-25
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
Gabriele Galli 2021-6-25
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

请先登录，再进行评论。

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

采纳的回答

Jan
Jan 2021-6-25
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.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by