专题3-4 插值与拟合+数值微分积分

1、专题3插值与拟合说明：专题 3作业插值与拟合-数学实验P81-实验内容 1（d）,3,5,6 1、p81-1-d（数值、画图都要体现出来）第一题d函数lagr文件function y=lagr1(x0,y0,x)n=length(x0);m=length(x);for i=1:m z=x(i); s=0; for k=1:n p=1; for j=1:n; if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end具体程序：x0=-2:2;%取n=5y0= exp(-x0.2);x=-2:0.1:2;%取m=41y

2、=exp(-x.2);y1=lagr1(x0,y0,x);y2=interp1(x0,y0,x);y3=spline(x0,y0,x);for k=1:11 xx(k)=x(19+2*k);%x>=0且间隔0.2的插值点 yy(k)=y(19+2*k);%x>=0且间隔0.2的函数值 yy1(k)=y1(19+2*k);%x>=0且间隔0.2的拉格朗日插值 yy2(k)=y2(19+2*k);%x>=0且间隔0.2的分段线性插值 yy3(k)=y3(19+2*k);%x>=0且间隔0.2的三次样条插值endxx;yy;yy1;yy2;yy3'z=0*x;p

3、lot(x,z,x,y,'b-', x,y1,'r')pauseplot(x,z,x,y,'k-', x,y2,'r')pauseplot(x,z,x,y,'k-',x,y3,'r')y=exp(-x2),-2<=x<=2,三种插值结果的计较ans = 0 1.0000 1.0000 1.0000 1.0000 0.2000 0.9608 0.9698 0.8736 0.9623 0.4000 0.8521 0.8815 0.7472 0.8617 0.6000 0.6977 0.742

4、7 0.6207 0.7168 0.8000 0.5273 0.5657 0.4943 0.5459 1.0000 0.3679 0.3679 0.3679 0.3679 1.2000 0.2369 0.1714 0.2980 0.2011 1.4000 0.1409 0.0036 0.2281 0.0642 1.6000 0.0773 -0.1035 0.1581 -0.0243 1.8000 0.0392 -0.1126 0.0882 -0.04572.0000 0.0183 0.0183 0.0183 0.0183一、题目：2000年世界人口预测60年代末世界人口增长情况：年： 1960

5、 1961 1962 1963 1964 1965 1966 1967 1968 人口：29.72 30.61 31.51 32.13 32.34 32.85 33.56 34.20 34.83记人口数为N(t)，采用指数函数N=exp(a+bt)对数据进行拟合，画出拟合曲线及散点图。(1)使用Matlab中的polyfit； (2) 直接使用最小二乘法求解（超定）方程组。2、 p81-3题程序：clc,clearx=1:0.5:10; y=x.3-6*x.2+5*x-3; y0=y+rand(1); figure(1)f1=polyfit(x,y0,3) %输出多项式系数y1=polyval

6、(f1,x); %计算各x点的拟合值plot(x,y,'+',x,y1) grid on title('三次拟合曲线'); c1=1 -6 5 -3'%原三次多项式系数xishubijiao3=c1,f1',c1-f1'%三次拟合系数比较wucha3=y',y0',y1',y1'-y',y1'-y0'%三次拟合误差figure(2); f2=polyfit(x,y0,2) %2次多项式拟合y2=polyval(f2,x); plot(x,y,'+',x,y2); gr

7、id on title('二次拟合曲线'); wucha2=y',y0',y2',y2'-y',y2'-y0'%二次拟合误差figure(3); f3=polyfit(x,y0,4) %4次多项式拟合y3=polyval(f3,x); plot(x,y,'+',x,y3) grid on title('四次拟合曲线');wucha4=y',y0',y3',y3'-y',y3'-y0'%四次拟合误差bjiao3_4=(wucha3(:,3

8、)'*(wucha3(:,3)-(wucha4(:,3)'*(wucha4(:,3)%比较三次四次的误差平方和（拟合数据和原来没有干扰的数据进行比较）结果：f1 = Columns 1 through 3 0.999999999999998 -5.999999999999961 4.999999999999779 Column 4 -2.450276391708465xishubijiao3 = 1.000000000000000 0.999999999999998 0.000000000000002 -6.000000000000000 -5.999999999999961

9、-0.000000000000039 5.000000000000000 4.999999999999779 0.000000000000221 -3.000000000000000 -2.450276391708465 -0.549723608291535wucha3 = 1.0e+002 * Columns 1 through 3 -0.030000000000000 -0.024502763917089 -0.024502763917087 -0.056250000000000 -0.050752763917089 -0.050752763917087 -0.09000000000000

10、0 -0.084502763917089 -0.084502763917088 -0.123750000000000 -0.118252763917089 -0.118252763917088 -0.150000000000000 -0.144502763917089 -0.144502763917088 -0.161250000000000 -0.155752763917089 -0.155752763917089 -0.150000000000000 -0.144502763917089 -0.144502763917089 -0.108750000000000 -0.1032527639

11、17089 -0.103252763917089 -0.030000000000000 -0.024502763917089 -0.024502763917089 0.093750000000000 0.099247236082911 0.099247236082911 0.270000000000000 0.275497236082911 0.275497236082911 0.506250000000000 0.511747236082911 0.511747236082911 0.810000000000000 0.815497236082911 0.815497236082911 1.

12、188750000000000 1.194247236082912 1.194247236082911 1.650000000000000 1.655497236082911 1.655497236082911 2.201250000000000 2.206747236082912 2.206747236082911 2.850000000000000 2.855497236082911 2.855497236082909 3.603750000000000 3.609247236082911 3.609247236082910 4.470000000000000 4.475497236082

13、911 4.475497236082909 Columns 4 through 5 0.005497236082913 0.000000000000002 0.005497236082913 0.000000000000001 0.005497236082912 0.000000000000001 0.005497236082912 0.000000000000000 0.005497236082912 0.000000000000000 0.005497236082911 0 0.005497236082911 -0.000000000000000 0.005497236082911 -0.

14、000000000000000 0.005497236082911 -0.000000000000000 0.005497236082911 -0.000000000000000 0.005497236082911 -0.000000000000000 0.005497236082911 -0.000000000000000 0.005497236082911 -0.000000000000000 0.005497236082911 -0.000000000000001 0.005497236082911 -0.000000000000001 0.005497236082911 -0.0000

15、00000000001 0.005497236082909 -0.000000000000002 0.005497236082911 -0.000000000000001 0.005497236082909 -0.000000000000002f2 = 10.499999999999993 -72.299999999999940 89.949723608291023wucha2 = 1.0e+002 * Columns 1 through 3 -0.030000000000000 -0.024502763917089 0.281497236082911 -0.056250000000000 -

16、0.050752763917089 0.051247236082911 -0.090000000000000 -0.084502763917089 -0.126502763917089 -0.123750000000000 -0.118252763917089 -0.251752763917089 -0.150000000000000 -0.144502763917089 -0.324502763917089 -0.161250000000000 -0.155752763917089 -0.344752763917089 -0.150000000000000 -0.14450276391708

17、9 -0.312502763917089 -0.108750000000000 -0.103252763917089 -0.227752763917088 -0.030000000000000 -0.024502763917089 -0.090502763917089 0.093750000000000 0.099247236082911 0.099247236082911 0.270000000000000 0.275497236082911 0.341497236082911 0.506250000000000 0.511747236082911 0.636247236082911 0.8

18、10000000000000 0.815497236082911 0.983497236082910 1.188750000000000 1.194247236082912 1.383247236082910 1.650000000000000 1.655497236082911 1.835497236082911 2.201250000000000 2.206747236082912 2.340247236082910 2.850000000000000 2.855497236082911 2.897497236082911 3.603750000000000 3.6092472360829

19、11 3.507247236082909 4.470000000000000 4.475497236082911 4.169497236082910 Columns 4 through 5 0.311497236082911 0.305999999999999 0.107497236082911 0.102000000000000 -0.036502763917089 -0.042000000000000 -0.128002763917089 -0.133500000000000 -0.174502763917089 -0.180000000000000 -0.183502763917089

20、-0.189000000000000 -0.162502763917089 -0.168000000000000 -0.119002763917088 -0.124500000000000 -0.060502763917089 -0.066000000000000 0.005497236082911 -0.000000000000000 0.071497236082911 0.066000000000000 0.129997236082911 0.124500000000000 0.173497236082910 0.167999999999999 0.194497236082910 0.18

21、8999999999999 0.185497236082911 0.179999999999999 0.138997236082910 0.133499999999999 0.047497236082910 0.041999999999999 -0.096502763917091 -0.102000000000002 -0.300502763917091 -0.306000000000002f3 = Columns 1 through 3 0.000000000000000 0.999999999999998 -5.999999999999990 Columns 4 through 5 4.9

22、99999999999997 -2.450276391708882wucha4 = 1.0e+002 * Columns 1 through 3 -0.030000000000000 -0.024502763917089 -0.024502763917089 -0.056250000000000 -0.050752763917089 -0.050752763917089 -0.090000000000000 -0.084502763917089 -0.084502763917089 -0.123750000000000 -0.118252763917089 -0.118252763917089

23、 -0.150000000000000 -0.144502763917089 -0.144502763917088 -0.161250000000000 -0.155752763917089 -0.155752763917088 -0.150000000000000 -0.144502763917089 -0.144502763917088 -0.108750000000000 -0.103252763917089 -0.103252763917088 -0.030000000000000 -0.024502763917089 -0.024502763917088 0.093750000000

24、000 0.099247236082911 0.099247236082912 0.270000000000000 0.275497236082911 0.275497236082912 0.506250000000000 0.511747236082911 0.511747236082912 0.810000000000000 0.815497236082911 0.815497236082911 1.188750000000000 1.194247236082912 1.194247236082912 1.650000000000000 1.655497236082911 1.655497

25、236082911 2.201250000000000 2.206747236082912 2.206747236082912 2.850000000000000 2.855497236082911 2.855497236082912 3.603750000000000 3.609247236082911 3.609247236082912 4.470000000000000 4.475497236082911 4.475497236082912 Columns 4 through 5 0.005497236082911 -0.000000000000000 0.005497236082911

26、 -0.000000000000000 0.005497236082911 -0.000000000000000 0.005497236082911 0.000000000000000 0.005497236082912 0.000000000000000 0.005497236082912 0.000000000000000 0.005497236082912 0.000000000000000 0.005497236082912 0.000000000000000 0.005497236082912 0.000000000000000 0.005497236082912 0.0000000

27、00000000 0.005497236082912 0.000000000000000 0.005497236082912 0.000000000000000 0.005497236082911 0 0.005497236082911 0 0.005497236082911 0 0.005497236082911 0 0.005497236082912 0.000000000000001 0.005497236082912 0.000000000000001 0.005497236082912 0.000000000000001bjiao3_4 = -5.238689482212067e-0

28、10图像：分析：在本题中利用了二次拟合，三次拟合，四次拟合分别拟合出函数y=x.3-6*x.2+5*x-3的多项式，并画出相应的拟合曲线。运行之后，对二次拟合曲线，三次拟合曲线，四次拟合曲线的图像进行比较分析，发现：相比较而言，三次拟合的吻合性更高，即精度更好。二次效果最差。3、p81 第五题：（要求程序、计算结果和图形）n=0000 1153 2045 2800 3466 4068 4621 5135 5619 6152;t=0 20 40 60 80 100 120 140 160 183.5;R=n'.2,n'aa=Rt'y=aa(1)*n'.2+aa(2

29、)*n'plot(n,t,'b+',n,y,'r')xlabel('n'),ylabel('T')t=0 20 40 60 80 100 120 140 160 184;n=0 1141 2019 2760 3413 4004 4545 5051 5525 6061;A=(n.2)' n'X= At' 注意：（1）、如果直接用polyfit做多项式拟合，会得出一般二次多项式；（2）、另外，可以做变换，两边同时除以n，然后用一次函数拟合。3、p81-6题（本答案在闾宇提供的基础上修改）分析：本题已经建

30、立好了模型！只需对其进行求解： 由题目所给的方程 （已知）可见，v(t)与（T）成指数变化关系，所以在通过曲线拟合的时候，使用指数曲线。（注：这是一种非线性曲线拟合，首先要进行变量代换，为方便计算做以下变量代换） 用代替，是拟合后的曲线方程 对 变形后取对数，有 。 令 则： 程序：t=0.5 1 2 3 4 5 7 9 ;v1=6.36 6.48 7.26 8.22 8.66 8.99 9.43 9.63;y=log(10-v1);a=polyfit(t,y,1);a1=a(1);a2=a(2);T=-1/a1v0=10-exp(a2)v2=10-(10-v0)*exp(-t/T);拟合后的

31、函数值plot(t,v1,'k+',t,v2,'r')结果：t0 = 3.5269v0 =5.6221专题4 数值积分微分专题4作业：数值积分与数值微分-数学实验P97-实验内容：2(d),3,9+补充完成p96人口增长率二、数值积分与微分部分0（未布置）97-1题利用三种公式计算积分，另外利用三点公式计算每点处的导数，数值计算结果都与精确值进行比较。(数值比较需列表、另外还要求画图比较)程序：syms xy=(x+sin(x/3);z=int(y,0.3,1.5)z=zxx=0.3:0.2:1.5;y=xx+sin(xx./3);z1=sum(y(1:6)*0

32、.2z2=sum(y(2:7)*0.2z3=quad('x+sin(x./3)',0.3,1.5)z4=trapz(xx,y)u1=z-z1,u2=z-z2,u3=z-z3,u4=z-z4h=0.2;x=0.3:h:1.5;y=0.3895 0.6598 0.9147 1.1611 1.3971 1.6212 1.8325;format longz1=sum(y(1:6)*h,z2=sum(y(2:7)*h,z3=trapz(y)*h,z4=quad('y',0.3,1.5),z5=quadl('y',0.3,1.5),format short1

33、、 p97-2题选择一些函数用梯形、辛普森和随机模拟三种方法计算积分。改变步长（对梯形公式），改变精度要求（对辛普森公式），改变随机点数目（对随机模拟），进行比较、分析。如下函数供选择参考：d.（1组同学提供答案）程序：编辑函数M文件：function y=fun(x)y=(1./sqrt(2*pi)*exp(-x.2./2)；程序：clc,clearformat longh=4/50; %步长x=-2:h:2;y=fun(x);z1=trapz(y)*h %梯形公式z2=quad('fun',-2,2) %辛普森公式n=1000; %随机模拟方法x1=rand(1,n);y1

34、=fun(x1.*2);z3=sum(y1)*4/n当对梯形公式取h=4/50,4/100,4/10000;对辛普森公式分别取精度为 ，对随机点数目取n=1000,10000,100000,得到的结果如下：梯形公式辛普森公式蒙特卡罗方法0.954384567677890.954499438241540.939992110595860.954470941689640.954499736107350.957943385948660.954499733224120.954491037397360.95427317381756分析：对于梯形公式，步长越小，结果越精确；对于辛普森公式，在一般情况精度下的

35、结果都很精确，辛普森有很好的优越性，但是需要函数表达式；蒙特卡罗方法，虽然结果具有随机性，但是随着n的增大，其计算结果越来越接近准确值。2.p98-3题（要求15位小数显示）x=0:0.1:1;y=(1-x.2).(0.5);z1=4*trapz(x,y)z2=4*quad('(1-x.2).(0.5)',0,1)z3=4*quad8('(1-x.2).(0.5)',0,1)n=10000000;x3=rand(1,n);yy=(1-x3.2).(0.5);z4=4*sum(yy)/n以下是闾宇同学1组（1）蒙特卡罗方法创建M文件function y=f(x)m

36、=0;for n=1:x if rand(1)2+rand(1)2<=1 m=m+1; end;end;4*m/x输入Format long>> f(10000);ans = 3.132400000000000>> f(50000);ans = 3.139920000000000>> f(100000);ans = 3.147320000000000（2）数值积分法 设 ，将区间0,1n等分，取 梯形法： 辛普森法： 创建M文件fun.mfunction y=fun(x)y=4./(1+x.2);创建M文件f(x)function y=f(x)for

37、n=1:x-1 a(n)=2*fun(n/x);end;vpa(1/(2*k)*(sum(a)+fun(0)+fun(1)输入：digits(30) %保留小数点后30位>> f(100)ans =3.14157598692312900467982217378>> f(500)ans =3.14159198692313035294887413329>> f(10000)ans =3.1415926519231400781961838223. p99-9题编程并最终回答问题。 题目：图7是一个国家的地图，为了算出它的国土面积，首先对地图作如下测量：以由西向东方

38、向为x轴，由南到北方向为y轴，选择方便的原点，并将从最西边界点到最东边界点在X轴上的区间适当地划分为若干段，在每个分点的y轴方向测出南边界点和北边界点的y坐标y1和y2，这样就得到了表中的数据（单位mm）。  根据地图的比例我们知道18mm相当于40km，试由测量数据计算该国土的近似面积，与它的精确值41288km2比较。（2012本-代雄，雷雨锦，李颖艺提供）程序1：>> x=7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 61.0 68.5 76.5 80.5 91.0 96.0 101.0 104.0 106.5 111.5

39、118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0;y1=44 45 47 50 50 38 30 30 34 36 34 41 45 46 43 37 33 28 32 65 55 54 52 50 66 66 68;y2=44 59 70 72 93 100 110 110 110 117 118 116 118 118 121 124 121 121 121 122 116 83 81 82 86 85 68;figure(1)plot(x,y1,'o')hold on plot(x,y2,'rp')hold o

40、nxi=10:10:100yi=interp1(x,y1,xi,'spline')yi=interp1(x,y2,xi,'spline')xx=min(x):0.1:max(x);yy1=interp1(x,y1,xx,'spline');plot(xx,yy1)hold onyy2=interp1(x,y2,xx,'spline');plot(xx,yy2,'r')hold onxi = 10 20 30 40 50 60 70 80 90 100yi = 44.6613 51.4824 53.2436 39.1

41、722 30.9541 35.8817 34.5265 44.6028 46.3324 38.2206yi = 56.3424 72.5187 87.6339 98.9170 108.8743 115.7368 117.3082 117.7779 117.8652 124.2328另外：两种方法，可以参考。结果：方法一本题利用蒙特卡罗方法计算国土面积，通过for循环产生随机点，if语句限定曲线所围区域内。 通过for循环多次实验，求的平均值。>> x=7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 61.0 68.5 76.5 80

42、.5 91.0 . 96.0 101.0 104.0 106.5 111.5 118.0 123.5 136.5 142.0 146.0 150.0 . 157.0 158.0; y1=44.0 45.0 47.0 50.0 50.0 38.0 30.0 30.0 34.0 36.0 34.0 41.0 45.0 46.0 . 43.0 37.0 33.0 28.0 32.0 65.0 55.0 54.0 52.0 50.0 66.0 66.0 68.0; y2=44.0 59.0 70.0 72.0 93.0 100.0 110.0 110.0 110.0 117.0 118.0 116.0

43、 . 118.0 118.0 121.0 124.0 121.0 121.0 121.0 122.0 116.0 83.0 81.0 . 82.0 86.0 85.0 68.0;L=max(x)-min(x);H=max(y2)-min(y1);s=0;for k=1:10 n=10000;m=0;for i=1:nu(i)=unifrnd(min(x),max(x);v(i)=unifrnd(min(y1),max(y2);f1=interp1(x,y1,u(i),'cubic');f2=interp1(x,y2,u(i),'cubic');if(v(i)&g

44、t;=f1&&v(i)<=f2)m=m+1;endendS(k)=1600*L*H*m/(n*182);s=s+S(k);end s=s/10s = 4.2552e+004方法二：x=7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 61.0 68.5 76.5 80.5 91.0 . 96.0 101.0 104.0 106.5 111.5 118.0 123.5 136.5 142.0 146.0 150.0 . 157.0 158.0; y1=44.0 45.0 47.0 50.0 50.0 38.0 30.0 30.0 34.

45、0 36.0 34.0 41.0 45.0 46.0 . 43.0 37.0 33.0 28.0 32.0 65.0 55.0 54.0 52.0 50.0 66.0 66.0 68.0; y2=44.0 59.0 70.0 72.0 93.0 100.0 110.0 110.0 110.0 117.0 118.0 116.0 . 118.0 118.0 121.0 124.0 121.0 121.0 121.0 122.0 116.0 83.0 81.0 . 82.0 86.0 85.0 68.0;newx=7:0.1:158;newy1=interp1(x,y1,newx,'lin

46、ear');newy2=interp1(x,y2,newx,'linear');area=sum(newy2-newy1)*0.1/182*1600)area = 4.2414e+004结果分析：对于方法一（插值更光滑），采用随机投点法结合分段三次插值模拟上下边界，求国土面积；方法二（简单，实用，不光滑），同样采取随机投点法，但结合分段线性插值模拟边界求解。可以看出方法二得出的结果要比方法一得出的结果相对精确值更接近些，说明国土边界是不规则，不光滑的。方法一、二都各有所长。补充4.完成p96人口增长率 要求：1）编程计算p96表格显示的增长率 2）编程计算1980年人口

47、x=1900:10:2000y=76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 281.4r(1)=(-3*y(1)+4*y(2)-y(3)/20/y(1)r(11)=(y(9)-4*y(10)+3*y(11)/20/y(11)for i=2:10 r(i)=(y(i+1)-y(i-1)/20/y(i)endr年份 1970 1972 1974 1976 1978 1980年增长率（%） 0.87 0.85 0.89 0.91 0.95 1.10为了从1.3表2估计1980年的人口，仍用上面的记号，人口增长率满足微分方程它在初

48、始条件下的解为 因为在题目中增长率以离散型数据给出（表3），所以上式要用数值积分计算.用梯形公式计算的结果是，1980年该地区人口为230.2（百万）.%p96 renkou(2)clc,cleary(1)=210r=0.87 0.85 0.89 0.91 0.95 1.10*0.01；t=0:2:10；for i=2:6 y(i)=y(1)*exp(2*trapz(r(1:i);endy'或者%p96 renkou(2)clc,cleary(1)=210r=0.87 0.85 0.89 0.91 0.95 1.10*0.01t=0:2:10for i=2:6 y(i)=y(1)*ex

49、p(trapz(t(1:i),r(1:i);endy'第七次 第四章Lingo练习程序1： max 72x1+64x2st2) x1+x2<503) 12x1+8x2<4804) 3x1<100Endmax 72x1+64x2st2) x1+x2<503) 12x1+8x2<4804) 3x1<100end LP OPTIMUM FOUND AT STEP 2 OBJECTIVE FUNCTION VALUE 1) 3360.000 VARIABLE VALUE REDUCED COST X1 20.000000 0.000000 X2 30.000000 0.000000 ROW SLACK OR SURPLUS DUAL PRICES 2) 0.000000 48.000000 3) 0.000000 2.0