I am confused about why my code doesn't lead to a fitting result.

1 次查看（过去 30 天）
Jinglei
Jinglei 2023-2-2
编辑： Walter Roberson 2023-2-6
Ran in:
x=[0 1 2 3 4 5 6 7 8 9 10]
x = 1×11
0 1 2 3 4 5 6 7 8 9 10
xs = 0:0.01:120;
y=[0 4 9 16 25 36 49 64 81 100 121]
y = 1×11
0 4 9 16 25 36 49 64 81 100 121
myfunc=@fitting;
rng('shuffle');
best_r = inf;
best_p = [-inf, -inf];
best_c = [-inf];
a=0.00001;
b=3;
d=2;
for iters = 1 : 100
p0 = a+(b-a).*rand(1,2);
c0 = a+(d-a).*rand(1,1);
ytotal0 = fitting(p0,c0,xs);
try
[p, c, r] =nlinfit(x,y,myfunc,p0,c0);
if any(p < 0), or (c < 0)
fprintf('rejected negative output on iteration %d\n', iters);
continue;
end
if any(norm(r)==best_r)
continue;
end
nr = norm(r);
if nr < best_r
fprintf('improved on iteration %d\n', iters);
best_r = nr;
best_p = p;
best_c = c;
end
catch ME
fprintf('failed iteration %d\n', iters);
end
end
failed iteration 1 failed iteration 2 failed iteration 3 failed iteration 4 failed iteration 5 failed iteration 6 failed iteration 7 failed iteration 8 failed iteration 9 failed iteration 10 failed iteration 11 failed iteration 12 failed iteration 13 failed iteration 14 failed iteration 15 failed iteration 16 failed iteration 17 failed iteration 18 failed iteration 19 failed iteration 20 failed iteration 21 failed iteration 22 failed iteration 23 failed iteration 24 failed iteration 25 failed iteration 26 failed iteration 27 failed iteration 28 failed iteration 29 failed iteration 30 failed iteration 31 failed iteration 32 failed iteration 33 failed iteration 34 failed iteration 35 failed iteration 36 failed iteration 37 failed iteration 38 failed iteration 39 failed iteration 40 failed iteration 41 failed iteration 42 failed iteration 43 failed iteration 44 failed iteration 45 failed iteration 46 failed iteration 47 failed iteration 48 failed iteration 49 failed iteration 50 failed iteration 51 failed iteration 52 failed iteration 53 failed iteration 54 failed iteration 55 failed iteration 56 failed iteration 57 failed iteration 58 failed iteration 59 failed iteration 60 failed iteration 61 failed iteration 62 failed iteration 63 failed iteration 64 failed iteration 65 failed iteration 66 failed iteration 67 failed iteration 68 failed iteration 69 failed iteration 70 failed iteration 71 failed iteration 72 failed iteration 73 failed iteration 74 failed iteration 75 failed iteration 76 failed iteration 77 failed iteration 78 failed iteration 79 failed iteration 80 failed iteration 81 failed iteration 82 failed iteration 83 failed iteration 84 failed iteration 85 failed iteration 86 failed iteration 87 failed iteration 88 failed iteration 89 failed iteration 90 failed iteration 91 failed iteration 92 failed iteration 93 failed iteration 94 failed iteration 95 failed iteration 96 failed iteration 97 failed iteration 98 failed iteration 99 failed iteration 100
c = best_c
c = -Inf
p = best_p
p = 1×2
-Inf -Inf
%c = lsqcurvefit(myfunc,c0,x,y)
ytotal = fitting(p,c,x);
Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t.
%for adding points
hold on
plot(x,y,'o', 'displayname', 'original')
hold off
xlim auto
ylim auto
legend show
%for adding points
hold on
xs = 0:0.1:150;
ytotal = fitting(p,c,xs);
Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t.
plot(xs,ytotal)
hold off
R2=1-(sum((fitting(p,c,x)-y).^2)/sum((y-mean(y)).^2));
Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t.
function ytotal = fitting(p,c,x)
xspan = x;
y0 = zeros(2,1);
[t,ye] = ode15s(@(x,y)eq2(x,y,p,c), xspan, y0);
%protect against failure of integration
if t(end) ~= xspan(end)
ytotal = inf(numel(xspan),1);
else
ytotal = ye(:,1)+ye(:,2)+c(1);
end
end
function dy=eq2(x,y,p,c)
%y(1)=CE,y(2)=LE
dy=zeros(2,1);
dy(1) = p(1) .* x;
dy(2) = p(2);
end

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

采纳的回答

Walter Roberson
Walter Roberson 2023-2-5
编辑：Walter Roberson 2023-2-5
Ran in:
I had to fix every different things to get it to work.
The improvement to summarize failed iterations instead of listing each one was cosmetic, not strictly required. It took a while before the iterations started working properly and I was tired of the long stream of fail messages.
x=[0 1 2 3 4 5 6 7 8 9 10]
x = 1×11
0 1 2 3 4 5 6 7 8 9 10
xs = 0:0.01:120;
y=[0 4 9 16 25 36 49 64 81 100 121]
y = 1×11
0 4 9 16 25 36 49 64 81 100 121
rng('shuffle');
best_r = inf;
best_p = [-inf, -inf];
best_c = [-inf];
a=0.00001;
b=3;
d=2;
arefailing = false;
firstfail = -inf;
currentfail = -inf;
for iters = 1 : 100
p0 = a+(b-a).*rand(1,2);
c0 = a+(d-a).*rand(1,1);
ytotal0 = fitting(p0,c0,xs);
try
[pc, r] = nlinfit(x,y(:),@(pc,x)fitting(pc(1:2),pc(3),x),[p0,c0]);
%with the try/catch we never reach this statement if the nlfit
%succeeded, so the fact we got here means that an iteration worked
%(but might have returned values we do not want.)
if arefailing
fprintf('failed iterations %d to %d\n', firstfail, currentfail);
arefailing = false;
end
p = pc(1:2); c = pc(3);
if any(p < 0), or (c < 0)
fprintf('rejected negative output on iteration %d\n', iters);
continue;
else
end
if any(norm(r)==best_r)
continue;
end
nr = norm(r);
if nr < best_r
fprintf('improved on iteration %d\n', iters);
best_r = nr;
best_p = p;
best_c = c;
end
catch ME
if arefailing
currentfail = iters;
else
arefailing = true;
firstfail = iters;
currentfail = iters;
disp(ME)
if ~isempty(ME.cause); disp(ME.cause{1}); end
end
end
end
improved on iteration 1 improved on iteration 4 improved on iteration 19
if arefailing
fprintf('failed iterations %d to %d\n', firstfail, currentfail);
if firstfail == 1 %everything failed
fprintf('failed all iterations!!! Not doing any plotting\n');
else
arefailing = false; %we got at least one valid iteration
end
end
if ~arefailing
c = best_c
p = best_p
%c = lsqcurvefit(myfunc,c0,x,y)
ytotal = fitting(p,c,x);
%for adding points
hold on
plot(x,y,'o', 'displayname', 'original')
hold off
xlim auto
ylim auto
legend show
%for adding points
hold on
xs = 0:0.1:150;
ytotal = fitting(p,c,xs);
plot(xs,ytotal)
hold off
R2=1-(sum((fitting(p,c,x)-y).^2)/sum((y-mean(y)).^2));
end
c = 0.4196
p = 1×2
1.9650 2.2203
function ytotal = fitting(p,c,x)
xspan = x;
y0 = zeros(2,1);
[t,ye] = ode15s(@(x,y)eq2(x,y,p,c), xspan, y0);
%protect against failure of integration
if t(end) ~= xspan(end)
ytotal = inf(numel(xspan),1);
else
ytotal = ye(:,1)+ye(:,2)+c;
end
end
function dy=eq2(x,y,p,c)
%y(1)=CE,y(2)=LE
dy=zeros(2,1);
dy(1) = p(1) .* x;
dy(2) = p(2);
end
  6 个评论
显示 4更早的评论
Walter Roberson
Walter Roberson 2023-2-6
In fact, it got stuck in the middle this time, running forever without stopping manually. Do you know what was wrong?
No. I do not know what is wrong with code that I have not seen.
Jinglei
Jinglei 2023-2-6
编辑：Walter Roberson 2023-2-6
The code can be seen as below.
x = [0
5
10
15
20
25
30
35
40
45
50
55
60
65
70
80
90
100
110
120
];
xs = 0:0.01:120;
y = [0
1.415699757
1.044160175
1.342848859
1.357419038
1.138866343
1.386559398
1.452125206
1.430269937
1.153436523
1.26271287
1.408414667
1.357419038
1.430269937
1.095155804
1.488550656
1.503120835
1.386559398
1.401129577
1.532261195
];
rng('shuffle');
best_r = inf;
best_p = [-inf, -inf,-inf,-inf,-inf,-inf,-inf,-inf];
best_c = [-inf,-inf];
a=0.00001;
b=1;
d=2;
arefailing = false;
firstfail = -inf;
currentfail = -inf;
for iters = 1 : 100
p0 = a+(b-a).*rand(1,8);
c0 = a+(d-a).*rand(1,2);
ytotal0 = fitting(p0,c0,xs);
try
[pc, r] = nlinfit(x,y(:),@(pc,x)fitting(pc(1:8),pc(9:10),x),[p0,c0]);
%with the try/catch we never reach this statement if the nlfit
%succeeded, so the fact we got here means that an iteration worked
%(but might have returned values we do not want.)
if arefailing
fprintf('failed iterations %d to %d\n', firstfail, currentfail);
arefailing = false;
end
p = pc(1:8); c = pc(9:10);
if any(p < 0), or (c < 0)
fprintf('rejected negative output on iteration %d\n', iters);
continue;
else
end
if any(norm(r)==best_r)
continue;
end
nr = norm(r);
if nr < best_r
fprintf('improved on iteration %d\n', iters);
best_r = nr;
best_p = p;
best_c = c;
end
catch ME
if arefailing
currentfail = iters;
else
arefailing = true;
firstfail = iters;
currentfail = iters;
disp(ME)
if ~isempty(ME.cause); disp(ME.cause{1}); end
end
end
end
if arefailing
fprintf('failed iterations %d to %d\n', firstfail, currentfail);
if firstfail == 1 %everything failed
fprintf('failed all iterations!!! Not doing any plotting\n');
else
arefailing = false; %we got at least one valid iteration
end
end
if ~arefailing
c = best_c
p = best_p
%c = lsqcurvefit(myfunc,c0,x,y)
ytotal = fitting(p,c,x);
%for adding points
hold on
plot(x,y,'o', 'displayname', 'original')
hold off
xlim auto
ylim auto
legend show
%for adding points
hold on
xs = 0:0.1:150;
ytotal = fitting(p,c,xs);
plot(xs,ytotal)
hold off
R2=1-(sum((fitting(p,c,x)-y).^2)/sum((y-mean(y)).^2));
end
function ytotal = fitting(p,c,x)
xspan = x;
y0 = zeros(4,1);
[t,ye] = ode15s(@(x,y)eq2(x,y,p,c), xspan, y0);
%protect against failure of integration
if t(end) ~= xspan(end)
ytotal = inf(numel(xspan),1);
else
ytotal = ye(:,1)+ye(:,2)+ye(:,3)+ye(:,4);
end
end
function dy=eq2(x,y,p,c)
%y(1)=CE,y(2)=LE
dy=zeros(4,1);
dy(1) = p(1) .* (2.5 - y(1)-y(2)-y(3)-y(4)) .* (c(1)-y(1)-y(2)-y(4)) - p(2) .* y(1);
dy(2) = p(3) .* y(1)-p(4).*y(2);
dy(3) = p(5).*(2.5-y(1)-y(2)-y(3)-y(4)).*(c(2)-y(3)-y(4))-p(6).*y(3);
dy(4) = p(7).*y(3).*(c(2)-y(1)-y(2)-y(4))-p(8).*y(4);
end

请先登录，再进行评论。

更多回答（1 个）

Bora Eryilmaz
Bora Eryilmaz 2023-2-2
If you print out the exception inside the catch block, it will tell you what the error is.
catch ME
ME
fprintf('failed iteration %d\n', iters);
To start with, you are not passing the right number of arguments to the fitting() function: it would be passed two arguments from nlinfit, but your code expects three input arguments.
  6 个评论
显示 4更早的评论
Jinglei
Jinglei 2023-2-5
Thank you for your suggestion. However, in this case, c and p are taken as un-defined parameters by Matlab. And I am confused how to fix that. Could you show me that?
Jinglei
Jinglei 2023-2-6
Thank you for your help and it does work this time. The only problem is the obtained parameters to some extent away from the optimal ones, which may not be a big problem for my data. Thank you for your patience.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Systems of Nonlinear Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by