fitting a curve using a complex equation as shown

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')

Answer 1

[Star Strider]

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')

.

[Vivek E K]

Thank you Strider. I'll try this algorithm.

[Star Strider]

As always, my pleasure!

That should work as it did here.

Use:

x0 = [-1000; 2000; -50; 1; 1];

to make the convergence more rapid.

The Global Optimization Toolbox also has solvers that will search the parameter space for solutions.

[Vivek E K]

Hello @Star Strider,

Happy New Year!

I was using the same cuve fit algorithm for a different set of data. But it is not giving a good fit. Can't we use the same algorithm for different data set of similar kind? or what may be the cause of it?

The data set used is given below.

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];

[Alex Sha]

It is more possible that the initial start-values of each parameter are inappropriate， of course, it is always a hard job for guessing proper initial start-values, refer to the result below:

Root of Mean Square Error (RMSE): 10386.9383983247
Sum of Squared Residual: 4423428060.91424
Correlation Coef. (R): 0.999943759470119
R-Square: 0.999887522103236

Parameter          	Best Estimate
--------------------	-------------
c1	-1007822382600.37
c2	-30889083.2743406
c3	1007827630712.89
c4	0.000518929889784698
c5	1.477652245583

Also, for your first data set, the result below is much better

Root of Mean Square Error (RMSE): 1.93524599034026
Sum of Squared Residual: 112.355311293842
Correlation Coef. (R): 0.997826661011103
R-Square: 0.995658045424567

Parameter          	Best Estimate
--------------------	-------------
c1	883.12162878574
c2	66.077201030121
c3	74.670493376168
c4	-0.0499468111882734
c5	0.383069709132556

[Star Strider]

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')

.