fitting a curve using a complex equation as shown

조회 수: 5(최근 30일)
Vivek E K
Vivek E K 2021년 6월 23일
댓글: Star Strider 2022년 1월 14일
Hi,
I am trying to fit a curve using this function by lsqnonlin.
But am not able to find a reliable fit using the code as shown below. Can anyone tell me how could i get a better fit.
Matlab Code
%%Function file
function F = f(c,w,stiff)
errors=c(1)+(c(2)*w.^2*c(4)^2./(1+w.^2*c(4)^2)+c(3)*w.^2*c(5)^2./(1+w.^2*c(5)^2))+1i*(c(2)*w*c(4)./(1+w.^2*c(4)^2)+c(3)*w*c(4)./(1+w.^2*c(5)^2))-stiff;
F = [real(errors(:)); imag(errors(:))];
%% Main program code
clc
clear all
format long
w=[1.02E+00,1.92E+00,2.94E+00,3.96E+00,4.98E+00,6.00E+00,7.01E+00,7.92E+00,8.94E+00,9.95E+00,1.10E+01,1.20E+01,1.30E+01,1.40E+01,1.49E+01,1.60E+01,1.70E+01,1.80E+01,1.90E+01,2.00E+01,2.10E+01,2.19E+01,2.30E+01,2.40E+01,2.50E+01,3.00E+01,3.50E+01,3.99E+01,4.50E+01,5.00E+01]';
stiff =[8.92E+02,9.12E+02,9.28E+02,9.37E+02,9.43E+02,9.50E+02,9.55E+02,9.59E+02,9.64E+02,9.68E+02,9.72E+02,9.74E+02,9.74E+02,9.78E+02,9.81E+02,9.83E+02,9.84E+02,9.87E+02,9.89E+02,9.91E+02,9.92E+02,9.94E+02,9.95E+02,9.94E+02,9.94E+02,1.00E+03,1.01E+03,1.01E+03,1.01E+03,1.02E+03]';
plot(w,stiff,'b*')
hold on
x0=[0 0 0 0 0]; % Arbitary inital guess
[x,resnorm] = lsqnonlin(@(c)f(c,w,stiff),x0);
c(1)=x(1);
c(2)=x(2);
c(3)=x(3);
c(4)=x(4);
c(5)=x(5);
z=c(1)+(c(2)*w.^2*c(4)^2./(1+w.^2*c(4)^2)+c(3)*w.^2*c(5)^2./(1+w.^2*c(5)^2))+1i*(c(2)*w*c(4)./(1+w.^2*c(4)^2)+c(3)*w*c(4)./(1+w.^2*c(5)^2));
z1=abs(z);
plot(w,z1,'r')

채택된 답변

Star Strider
Star Strider 2021년 6월 23일
Use lsqcurvefit instead. It does all the ‘heavy lifting’.
Try this —
format long G
f = @(c,w) c(1)+(c(2)*w.^2*c(4)^2./(1+w.^2*c(4)^2)+c(3)*w.^2*c(5)^2./(1+w.^2*c(5)^2))+1i*(c(2)*w*c(4)./(1+w.^2*c(4)^2)+c(3)*w*c(4)./(1+w.^2*c(5)^2));
w=[1.02E+00,1.92E+00,2.94E+00,3.96E+00,4.98E+00,6.00E+00,7.01E+00,7.92E+00,8.94E+00,9.95E+00,1.10E+01,1.20E+01,1.30E+01,1.40E+01,1.49E+01,1.60E+01,1.70E+01,1.80E+01,1.90E+01,2.00E+01,2.10E+01,2.19E+01,2.30E+01,2.40E+01,2.50E+01,3.00E+01,3.50E+01,3.99E+01,4.50E+01,5.00E+01]';
stiff =[8.92E+02,9.12E+02,9.28E+02,9.37E+02,9.43E+02,9.50E+02,9.55E+02,9.59E+02,9.64E+02,9.68E+02,9.72E+02,9.74E+02,9.74E+02,9.78E+02,9.81E+02,9.83E+02,9.84E+02,9.87E+02,9.89E+02,9.91E+02,9.92E+02,9.94E+02,9.95E+02,9.94E+02,9.94E+02,1.00E+03,1.01E+03,1.01E+03,1.01E+03,1.02E+03]';
x0 = rand(5,1);
[x,resnorm] = lsqcurvefit(@(c,w)abs(f(c,w)), x0, w, stiff)
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.
x = 5×1
-887.551435061344 1927.38724482353 -38.0486415450175 0.305270145980733 -0.0245243104477786
resnorm =
153.990966954166
figure
plot(w, stiff, 'pb')
hold on
plot(w, abs(f(x,w)), '-r')
hold off
grid
xlabel('w')
ylabel('stiff')
.
  댓글 수: 5
Star Strider
Star Strider 2022년 1월 14일
One problem I immediately noticed is that in the previous problelm, ‘w’ began at 0 whiule here it begins at 300. I have no idea if that makes any difference, however I created a second funciton that incorpirates an ‘offset’ for ‘w’ that can be estimated as an additional parameter.
However, I cannot get any sort of good fit on these data even using Global Optimization Toolbox functions. (I tried GlobalSearch and MultiStart first, then returned to ga since I have more experience with it.)
I will come back to this later to see if I can get anywhere with it.
The code I am using —
w=[300 305 310 315 320 325 330 335 340 345 350 355 360 365 370 375 380 385 390 395 400 405 410 415 420 425 430 435 440 445 450 455 460 465 470 475 480 485 490 495 500];
stiff =[3932520.873 4015693.586 4097988.694 4181314.464 4253465.759 4328014.218 4406125.004 4484569.939 4567617.217 4660316.661 4756506.294 4858930.943 4960939.872 5059458.756 5147133.194 5233133.891 5321741.762 5406976.925 5493062.328 5578653.654 5656765.584 5738375.124 5821908.655 5898492.521 5977227.204 6057701.956 6139178.58 6219528.066 6297853.615 6376801.961 6455759.104 6528528.198 6598112.246 6675012.222 6743401.917 6820245.288 6894851.756 6969871.718 7046138.351 7116057.231 7190824.713];
f = @(c,w) c(1)+(c(2)*w.^2*c(4)^2./(1+w.^2*c(4)^2)+c(3)*w.^2*c(5)^2./(1+w.^2*c(5)^2))+1i*(c(2)*w*c(4)./(1+w.^2*c(4)^2)+c(3)*w*c(4)./(1+w.^2*c(5)^2));
fws = @(c,w) c(1)+(c(2)*(w-c(6)).^2*c(4)^2./(1+(w-c(6)).^2*c(4)^2)+c(3)*(w-c(6)).^2*c(5)^2./(1+(w-c(6)).^2*c(5)^2))+1i*(c(2)*(w-c(6))*c(4)./(1+(w-c(6)).^2*c(4)^2)+c(3)*(w-c(6))*c(4)./(1+(w-c(6)).^2*c(5)^2));
ftns = @(c) norm(stiff - fws(c,w));
% x0 = [-1E3; 2E3; -50; 1; 1];
format longG
PopSz = 50;
Parms = 6;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms).*[randn*1E+8 randn*1E+4 randn*1E-2 randn*1E-3 randn*1E-3 1E-1], 'MaxGenerations',2E4, 'FunctionTolerance',1E-8);%, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
[x, fv] = ga(ftns, Parms, [],[],[],[],[-Inf(Parms-1,1); 0],[Inf(Parms-1,1); 500],[],[],opts)
Optimization terminated: maximum number of generations exceeded.
x = 1×6
1.0e+00 * 6635.48597835566 6618.92934084008 6632.06049319454 2413.97605039813 2377.32223536104 499.59111192986
fv =
36332739.0677012
Paramv = x(:)
Paramv = 6×1
1.0e+00 * 6635.48597835566 6618.92934084008 6632.06049319454 2413.97605039813 2377.32223536104 499.59111192986
figure
plot(w, stiff, 'pb')
hold on
plot(w, abs(fws(x,w)), '-r')
hold off
grid
xlabel('w')
ylabel('stiff')
.

댓글을 달려면 로그인하십시오.

추가 답변(0개)

Community Treasure Hunt

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

Start Hunting!

Translated by