k7=[0.5 30];
m7=lsqcurvefit(f,k7,t,M7) z7=f(m7,t) k8=[0.5 30];
m8=lsqcurvefit(f,k8,t,M8) z8=f(m8,t) k9=[0.5 30];
m9=lsqcurvefit(f,k9,t,M9) z9=f(m9,t)
plot(t2,M1,'b.'),hold
on,plot(t2,z1,'r-'),plot(t5,M2,'b.'),hold
on,plot(t5,z2,'r--'),plot(t3,M3,'b.'),hold on,plot(t3,z3,'y-'),plot(t4,M4,'b.'),hold on,plot(t4,z4,'y--'),plot(t,M5,'b.'),hold on,plot(t,M6,'b.'),hold
on,plot(t,z5,'g-'),hold
on,plot(t,z6,'g--'),plot(t,M7,'b.'),hold
on,plot(t,z8,'m--'),hold
on,plot(t,z7,'m-'),plot(t,M8,'b.'),hold on,plot(t,M9,'b.'),hold on,plot(t,z9,'k-') title('M取不同的值时每年增量的分布情况')
其运行结果如下：
Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun m1 =
46.0000 617.6667 z1 =
663.6667 709.6667 755.6667 Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun m2 =
46.6000 548.5000 z2 =
595.1000 641.7000 688.3000 734.9000 Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun m3 =
42.6857 477.2667 z3 =
519.9524 562.6381 605.3238 648.0095 690.6952 733.3810 Optimization terminated successfully:
Norm of the current step is less than OPTIONS.TolX m4 =
38.8810 408.5357 z4 =
Columns 1 through 7
447.4167 486.2976 525.1786 564.0595 602.9405 641.8214 680.7024 Column 8 719.5833
Optimization terminated successfully:
Norm of the current step is less than OPTIONS.TolX m5 =
35.6424 337.8667 z5 =
Columns 1 through 7
373.5091 409.1515 444.7939 480.4364 516.0788 551.7212 587.3636 Columns 8 through 10 623.0061 658.6485 694.2909 Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun m6 =
28.6121 281.1333 z6 =
Columns 1 through 7
309.7455 338.3576 366.9697 395.5818 424.1939 452.8061 481.4182 Columns 8 through 10 510.0303 538.6424 567.2545 Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun m7 =
22.8061 218.0667 z7 =
Columns 1 through 7
240.8727 263.6788 286.4848 309.2909 332.0970 354.9030 377.7091 Columns 8 through 10 400.5152 423.3212 446.1273 Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun m8 =
16.4364 156.4000 z8 =
Columns 1 through 7
172.8364 189.2727 205.7091 222.1455 238.5818 255.0182 271.4545 Columns 8 through 10 287.8909 304.3273 320.7636 Optimization terminated successfully:
First-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected m9 =
10.7636 100.2000 z9 =
Columns 1 through 7
110.9636 121.7273 132.4909 143.2545 154.0182 164.7818 175.5455 Columns 8 through 10 186.3091 197.0727 207.8364
附录11（M、D关系的稳定性检验）
x=[2000 3000 4000 5000 6000 7000 8000 9000 10000];
y=[10.7636 16.4364 22.8061 28.6121 35.6424 37.5357 40.4000 46.6000 46.0000]; f=inline('k(1)*x+k(2)','k','x'); k=[0.5 30];
m=lsqcurvefit(f,k,x,y) z=f(m,x)
plot(x,y,'.'),hold on,plot(x,z,'-') 其输出结果如下：
>> Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun m =
0.0046 4.0893 z =
Columns 1 through 7
13.2742 17.8666 22.4591 27.0516 31.6440 36.2365 40.8290 Columns 8 through 9 45.4214 50.0139
附录12（问题四稳定性分析）
>> x=[2000 3000 4000 5000 6000 7000 8000 9000 10000]'; >> X=[ones(9,1)x];
Warning: Future versions of MATLAB will require that whitespace, a comma, or a semicolon separate elements of a matrix. Please type
\>> Y=[10.7636 16.4364 22.8061 28.6121 35.6424 38.8810 42.6857 46.6000 46.0000]';
>> [b,bint,r,rint,stats]=regress(Y,X); >> b,bint,stats b = 3.9010 0.0047 bint =
-1.5746 9.3766 0.0039 0.0055 stats =
0.9616 175.1021 0.0000 >> rcoplot(r,rint) >> z=b(1)+b(2)*x z = 13.2832 17.9743 22.6653 27.3564 32.0475 36.7386 41.4296 46.1207 50.8118
>> plot(x,Y,'k+',x,z,'r')