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

专题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=5 y0= exp(-x0.^2);

x=-2:0.1:2;%取m=41 y=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的三次样条插值 end

[xx;yy;yy1;yy2;yy3]' z=0*x;

plot(x,z,x,y,'b--', x,y1,'r') pause

plot(x,z,x,y,'k--', x,y2,'r') pause

plot(x,z,x,y,'k--',x,y3,'r')

y=exp(-x^2),-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.7427 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.0457

2.0000 0.0183 0.0183 0.0183 0.0183

一、题目：2000年世界人口预测 60年代末世界人口增长情况：

年： 1960 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,clear x=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(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); grid 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))'*(wucha3(:,3))-(wucha4(:,3))'*(wucha4(:,3))%比较三次四次的误差[平方和]（拟合数据和原来没有干扰的数据进行比较）

结果： f1 =

Columns 1 through 3

0.999999999999998 -5.999999999999961 4.999999999999779

Column 4

-2.450276391708465

xishubijiao3 =

1.000000000000000 0.999999999999998 0.000000000000002 -6.000000000000000 -5.999999999999961 -0.000000000000039 5.000000000000000 4.999999999999779 0.000000000000221 -3.000000000000000 -2.450276391708465 -0.549723608291535

wucha3 =

1.0e+002 *

Columns 1 through 3

-0.030000000000000 -0.024502763917089 -0.024502763917087 -0.056250000000000 -0.050752763917089 -0.050752763917087 -0.090000000000000 -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.103252763917089 -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.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.475497236082911 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.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.000000000000001 0.005497236082909 -0.000000000000002 0.005497236082911 -0.000000000000001 0.005497236082909 -0.000000000000002

f2 =

10.499999999999993 -72.299999999999940 89.949723608291023

wucha2 =

1.0e+002 *

Columns 1 through 3

-0.030000000000000 -0.024502763917089 0.281497236082911 -0.056250000000000 -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.144502763917089 -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.810000000000000 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.609247236082911 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 -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.188999999999999 0.185497236082911 0.179999999999999 0.138997236082910 0.133499999999999 0.047497236082910 0.041999999999999 -0.096502763917091 -0.102000000000002 -0.300502763917091 -0.306000000000002 f3 =

Columns 1 through 3

0.000000000000000 0.999999999999998 -5.999999999999990

Columns 4 through 5

4.999999999999997 -2.450276391708882

wucha4 =

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 -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.093750000000000 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.655497236082911 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 -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.000000000000000 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.000000000000001 bjiao3_4 =

-5.238689482212067e-010

图像：

三次拟合曲线450400350300250200150100500-5012345678910

二次拟合曲线450400350300250200150100500-5012345678910

四次拟合曲线450400350300250200150100500-5012345678910

分析：在本题中利用了二次拟合，三次拟合，四次拟合分别拟合出函数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=R\\t'

y=aa(1)*n'.^2+aa(2)*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= A\\t'

注意： （1）、如果直接用polyfit做多项式拟合，会得出一般二次多项式； （2）、另外，可以做变换，两边同时除以n，然后用一次函数拟合。 3、p81-6题

（本答案在闾宇提供的基础上修改）

分析：本题已经建立好了模型！只需对其进行求解：

由题目所给的方程v(t)?v?(v?v0)exp(?t/?) （已知v?10）可见，v(t)与?（T）成指数变化关系，所以在通过曲线拟合的时候，使用指数曲线y?a2ea1x。（注：这是一种非线性曲线拟合，首先要进行变量代换，为方便计算做以下变量代换） ①用v1代替v(t)，v2是拟合后的曲线方程 ②对v(t)?10?(10?v0)exp(?t/?) 变形后取对数，有

ln(10?v(t))?ln(10?v0)?(?t/?) 。

③ 令y?ln(10?v(t)),a2?ln(10?v0),a1??1/?

则：v0?10?exp(a2),???1/a1

程序：

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/a1

v0=10-exp(a2)

v2=10-(10-v0)*exp(-t/T);拟合后的函数值 plot(t,v1,'k+',t,v2,'r')

结果： t0 =

3.5269 v0 =

5.6221

109.598.587.576.560123456789

专题4 数值积分微分

专题4作业：数值积分与数值微分-《数学实验》P97-实验内容：2(d),3,9+补充完成p96人口增长率

二、数值积分与微分部分

0（未布置）97-1题利用三种公式计算积分，另外利用三点公式计算每点处的导数，数值计算结果都与精确值进行比较。(数值比较需列表、另外还要求画图比较)

程序： syms x

y=(x+sin(x/3)); z=int(y,0.3,1.5) z=z

xx=0.3:0.2:1.5; y=xx+sin(xx./3); z1=sum(y(1:6))*0.2 z2=sum(y(2:7))*0.2

z3=quad('x+sin(x./3)',0.3,1.5) z4=trapz(xx,y)

u1=z-z1,u2=z-z2,u3=z-z3,u4=z-z4

h=0.2;x=0.3:h:1.5;

y=[0.3895 0.6598 0.9147 1.1611 1.3971 1.6212 1.8325]; format long

z1=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 short

1、p97-2题

选择一些函数用梯形、辛普森和随机模拟三种方法计算积分。改变步长（对梯形公式），改变精度要求（对辛普森公式），改变随机点数目（对随机模拟），进行比较、分析。如下函数供选择参考：

1?x2d.y?e,?2?x?2.

2?（1组同学提供答案） 程序：编辑函数M文件： function y=fun(x)

y=(1./sqrt(2*pi))*exp(-x.^2./2)；

程序： clc,clear

format long h=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=fun(x1.*2); z3=sum(y1)*4/n

当对梯形公式取h=4/50,4/100,4/10000;对辛普森公式分别取精度为

?310,1?708 1，对随机点数目取n=1000,10000,100000,得到的结果如下： ,?02梯形公式 辛普森公式 蒙特卡罗方法 0.95438456767789 0.95449943824154 0.93999211059586 0.95447094168964 0.95449973610735 0.95794338594866 0.95449973322412 0.95449103739736 0.95427317381756 分析：对于梯形公式，步长越小，结果越精确；对于辛普森公式，在一般情况精度下的结果都很精确，辛普森有很好的优越性，但是需要函数表达式；蒙特卡罗方法，虽然结果具有随机性，但是随着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=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）数值积分法

A?4?1dx?? 01?x21设y(x)?1k1 ，将区间[0,1]n等分，取 x?,y?kk21?x2n1?xk2梯形法：A?[2(y1?y2?????yn?1)?y0?yn]

n辛普森法：

1[(y0?y2m)?2(y2?y4?????y2m?2)?4(y1?y3?????y2m?1)] 6m创建M文件fun.m function y=fun(x) y=4./(1+x.^2); 创建M文件f(x) function y=f(x) for 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.141592651923140078196183822

3.p99-9题编程并最终回答问题。

题目：图7是一个国家的地图，为了算出它的国土面积，首先对地图作如下测量：以由西向东方向为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 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 on

xi=10:10:100

yi=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 on

yy2=interp1(x,y2,xx,'spline'); plot(xx,yy2,'r') hold on xi =

10 20 30 40 50 60 70 80 90 100 yi =

44.6613 51.4824 53.2436 39.1722 30.9541 35.8817 34.5265 46.3324 38.2206 yi =

56.3424 72.5187 87.6339 98.9170 108.8743 115.7368 117.3082 117.8652 124.2328

14012010080604020020406080100120140160

44.6028 117.7779 另外：两种方法，可以参考。

结果： 方法一

本题利用蒙特卡罗方法计算国土面积，通过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.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 ... 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:n

u(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)>=f1&&v(i)<=f2) m=m+1; end end

S(k)=1600*L*H*m/(n*18^2); s=s+S(k); end s=s/10 s =

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.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,'linear'); newy2=interp1(x,y2,newx,'linear');

area=sum((newy2-newy1)*0.1/18^2*1600)

area =

4.2414e+004

结果分析：对于方法一（插值更光滑），采用随机投点法结合分段三次插值模拟上下边界，求国土面积；方法二（简单，实用，不光滑），同样采取随机投点法，但结合分段线性插值模拟边界求解。可以看出方法二得出的结果要比方法一得出的结果相对精确值更接近些，说明国土边界是不规则，不光滑的。方法一、二都各有所长。

补充4.完成p96人口增长率

要求：1）编程计算p96表格显示的增长率 2）编程计算1980年人口

x=1900:10:2000

y=[76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 281.4]

r(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) end r

年份 1970 1972 1974 1976 1978 1980 年增长率（%） 0.87 0.85 0.89 0.91 0.95 1.10 为了从1.3表2估计1980年的人口，仍用上面的记号，人口增长率满足微分方程

dx?r(t)x(t),它在初始条件x(0)?x0下的解为 dtt0 x(t)?x0e?r(u)du

因为在题目中增长率r(u)以离散型数据给出（表3），所以上式要用数值积分计算.用梯形公式计算的结果是，1980年该地区人口为230.2（百万）.

%p96 renkou(2) clc,clear y(1)=210

r=[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))); end y' 或者

%p96 renkou(2) clc,clear

y(1)=210

r=[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(trapz(t(1:i),r(1:i))); end y'

第七次 第四章Lingo练习

程序1： max 72x1+64x2 st

2) x1+x2<50

3) 12x1+8x2<480 4) 3x1<100 End

max 72x1+64x2 st

2) x1+x2<50

3) 12x1+8x2<480 4) 3x1<100 end

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.000000 4) 40.000000 0.000000

NO. ITERATIONS= 2

RANGES IN WHICH THE BASIS IS UNCHANGED:

OBJ COEFFICIENT RANGES

VARIABLE CURRENT ALLOWABLE ALLOWABLE COEF INCREASE DECREASE X1 72.000000 24.000000 8.000000 X2 64.000000 8.000000 16.000000

RIGHTHAND SIDE RANGES

ROW CURRENT ALLOWABLE ALLOWABLE RHS INCREASE DECREASE 2 50.000000 10.000000 6.666667 3 480.000000 53.333332 80.000000 4 100.000000 INFINITY 40.000000 程序2：max 24x1+16x2+44x3+32x4-3x5-3x6 s.t.

2) 4x1+3x2+4x5+3x6<600 3) 4x1+2x2+6x5+4x6<480 4) x1+x5<100 5) x3-0.8x5=0 6) x4-0.75x6=0 End

LP OPTIMUM FOUND AT STEP 2

OBJECTIVE FUNCTION VALUE

1) 3460.800

VARIABLE VALUE REDUCED COST X1 0.000000 1.680000 X2 168.000000 0.000000 X3 19.200001 0.000000 X4 0.000000 0.000000 X5 24.000000 0.000000 X6 0.000000 1.520000

ROW SLACK OR SURPLUS DUAL PRICES 2) 0.000000 3.160000 3) 0.000000 3.260000 4) 76.000000 0.000000 5) 0.000000 44.000000 6) 0.000000 32.000000

NO. ITERATIONS= 2

RANGES IN WHICH THE BASIS IS UNCHANGED:

OBJ COEFFICIENT RANGES

VARIABLE CURRENT ALLOWABLE ALLOWABLE COEF INCREASE DECREASE X1 24.000000 1.680000 INFINITY X2 16.000000 8.150000 2.100000 X3 44.000000 19.750002 3.166667 X4 32.000000 2.026667 INFINITY X5 -3.000000 15.800000 2.533334 X6 -3.000000 1.520000 INFINITY

RIGHTHAND SIDE RANGES

ROW CURRENT ALLOWABLE ALLOWABLE RHS INCREASE DECREASE 2 600.000000 120.000000 280.000000 3 480.000000 253.333328 80.000000 4 100.000000 INFINITY 76.000000 5 0.000000 INFINITY 19.200001 6 0.000000 INFINITY 0.000000 程序3：max 2x1+3x2+4x3 s.t.

2) 1.5x1+3x2+5x3<600

3) 280x1+250x2+400x3<60000 End

LP OPTIMUM FOUND AT STEP 1

OBJECTIVE FUNCTION VALUE

1) 632.2581

VARIABLE VALUE REDUCED COST X1 64.516129 0.000000 X2 167.741928 0.000000 X3 0.000000 0.946237

ROW SLACK OR SURPLUS DUAL PRICES 2) 0.000000 0.731183 3) 0.000000 0.003226

NO. ITERATIONS= 1

RANGES IN WHICH THE BASIS IS UNCHANGED:

OBJ COEFFICIENT RANGES

VARIABLE CURRENT ALLOWABLE ALLOWABLE COEF INCREASE DECREASE X1 2.000000 1.360000 0.500000 X2 3.000000 1.000000 0.550000 X3 4.000000 0.946237 INFINITY

RIGHTHAND SIDE RANGES

ROW CURRENT ALLOWABLE ALLOWABLE RHS INCREASE DECREASE 2 600.000000 119.999992 278.571411 3 60000.000000 52000.000000 10000.000000 程序4：max 2x1+3x2+4x3 st

2) 1.5x1+3x2+5x3<600

3) 280x1+250x2+400x3<60000 end gin3

LP OPTIMUM FOUND AT STEP 0

OBJECTIVE FUNCTION VALUE

1) 632.2581

VARIABLE VALUE REDUCED COST X1 64.516129 0.000000 X2 167.741928 0.000000 X3 0.000000 0.946237

ROW SLACK OR SURPLUS DUAL PRICES 2) 0.000000 0.731183 3) 0.000000 0.003226

NO. ITERATIONS= 0

RANGES IN WHICH THE BASIS IS UNCHANGED:

OBJ COEFFICIENT RANGES

VARIABLE CURRENT ALLOWABLE ALLOWABLE COEF INCREASE DECREASE X1 2.000000 1.360000 0.500000 X2 3.000000 1.000000 0.550000 X3 4.000000 0.946237 INFINITY

RIGHTHAND SIDE RANGES

ROW CURRENT ALLOWABLE ALLOWABLE RHS INCREASE DECREASE 2 600.000000 119.999992 278.571411 3 60000.000000 52000.000000 10000.000000 程序5：model:

Max=4.8*x11+4.8*x21+5.6*x12+5.6*x22-10*x1-8*x2-6*x3; x11+x120; 0.4*x12-0.6*x22>0; x=x1+x2+x3; (x1-500)*x2=0; (x2-500)*x3=0; x1<500; x2<500; x3<500;

x>0; x11>0; x12>0; x21>0; x22>0; x1>0; x2>0; x3>0; end

Rows= 19 Vars= 8 No. integer vars= 0

Nonlinear rows= 2 Nonlinear vars= 3 Nonlinear constraints= 2

Nonzeros= 42 Constraint nonz= 28 Density=0.246

No. < : 5 No. =: 3 No. > : 10, Obj=MAX Single cols= 0

Optimal solution found at step: 4 Objective value: 4800.000

Variable Value Reduced Cost X11 500.0000 0.0000000 X21 500.0000 0.0000000 X12 0.0000000 0.0000000 X22 0.0000000 0.0000000 X1 0.1021405E-13 10.00000 X2 0.0000000 8.000000 X3 0.0000000 6.000000 X 0.0000000 0.0000000

Row Slack or Surplus Dual Price 1 4800.000 1.000000 2 0.0000000 9.600000 3 500.0000 0.0000000 4 0.0000000 -9.600000 5 0.0000000 -9.333333 6 0.0000000 0.0000000 7 0.0000000 0.0000000 8 0.0000000 0.0000000 9 500.0000 0.0000000 10 500.0000 0.0000000 11 500.0000 0.0000000 12 0.0000000 9.600000 13 500.0000 0.0000000

14 0.0000000 -0.2666667 15 500.0000 0.0000000 16 0.0000000 0.0000000 17 0.0000000 0.0000000 18 0.0000000 0.0000000 19 0.0000000 0.0000000

程序6：max 4.8x11+4.8x21+5.6x12+5.6x22-10x1-8x2-6x3 st

x-x1-x2-x3=0 x11+x12-x<500 x21+x22<1000 0.5x11-0.5x21>0 0.4x12-0.6x22>0 x1-500y1<0 x2-500y2<0 x3-500y3<0 x1-500y2>0 x2-500y3>0 end int y1 int y2 int y3

ENUMERATION COMPLETE. BRANCHES= 2 PIVOTS=

LAST INTEGER SOLUTION IS THE BEST FOUND RE-INSTALLING BEST SOLUTION...

OBJECTIVE FUNCTION VALUE

1) 5000.000

VARIABLE VALUE REDUCED COST Y1 1.000000 0.000000 Y2 1.000000 2000.000000 Y3 1.000000 1000.000000 X11 0.000000 1.200000 X21 0.000000 0.200000 X12 1500.000000 0.000000 X22 1000.000000 0.000000 X1 500.000000 0.000000 X2 500.000000 0.000000 X3 0.000000 0.000000

28 X 1000.000000 0.000000

ROW SLACK OR SURPLUS DUAL PRICES 2) 0.000000 6.000000 3) 0.000000 6.000000 4) 0.000000 5.000000 5) 0.000000 0.000000 6) 0.000000 -1.000000 7) 0.000000 0.000000 8) 0.000000 0.000000 9) 500.000000 0.000000 10) 0.000000 -4.000000 11) 0.000000 -2.000000

NO. ITERATIONS= 35

BRANCHES= 2 DETERM.= 1.000E 0

程序7：

Min 3x1+x2+3x3+3x4+x5+x6+3x7 s.t.

4x1+3x2+2x3+x4+x5>=50 x2+2x4+x5+3x6>=20 x3+x5+2x7>=15 end gin7

LP OPTIMUM FOUND AT STEP 4

OBJECTIVE FUNCTION VALUE

1) 26.66667

VARIABLE VALUE REDUCED COST X1 0.000000 1.666667 X2 11.666667 0.000000 X3 0.000000 1.666667 X4 0.000000 2.666667 X5 15.000000 0.000000 X6 0.000000 1.000000 X7 0.000000 1.666667

ROW SLACK OR SURPLUS DUAL PRICES 2) 0.000000 -0.333333 3) 6.666667 0.000000

4) 0.000000 -0.666667

NO. ITERATIONS= 4

RANGES IN WHICH THE BASIS IS UNCHANGED:

OBJ COEFFICIENT RANGES

VARIABLE CURRENT ALLOWABLE ALLOWABLE COEF INCREASE DECREASE X1 3.000000 INFINITY 1.666667 X2 1.000000 1.250000 1.000000 X3 3.000000 INFINITY 1.666667 X4 3.000000 INFINITY 2.666667 X5 1.000000 0.833333 0.666667 X6 1.000000 INFINITY 1.000000 X7 3.000000 INFINITY 1.666667

RIGHTHAND SIDE RANGES

ROW CURRENT ALLOWABLE ALLOWABLE RHS INCREASE DECREASE 2 50.000000 INFINITY 19.999998 3 20.000000 6.666667 INFINITY 4 15.000000 35.000000 9.999999

Model:

min=x1+x2+x3;

x1*r11+x2*r12+x3*r13>=50; x1*r21+x2*r22+x3*r23>=10; x1*r31+x2*r32+x3*r33>=20; x1*r41+x2*r42+x3*r43>=15; 4*r11+5*r21+6*r31+8*r41<=19; 4*r12+5*r22+6*r32+8*r42<=19; 4*r13+5*r23+6*r33+8*r43<=19; 4*r11+5*r21+6*r31+8*r41>=16; 4*r12+5*r22+6*r32+8*r42>=16; 4*r13+5*r23+6*r33+8*r43>=16; x1+x2+x3>=26; x1+x2+x3<=31; x1>=x2; x2>=x3;

@gin(x1);@gin(x2);@gin(x3); @gin(r11);@gin(r12);@gin(r13); @gin(r21);@gin(r22);@gin(r23); @gin(r31);@gin(r32);@gin(r33);

@gin(r41);@gin(r42);@gin(r43); end

Rows= 15 Vars= 15 No. integer vars= 15

Nonlinear rows= 4 Nonlinear vars= 15 Nonlinear constraints= 4

Nonzeros= 73 Constraint nonz= 58 Density=0.304

No. < : 4 No. =: 0 No. > : 10, Obj=MIN Single cols= 0

Optimal solution found at step: 28970 Objective value: 28.00000 Branch count: 2946

Variable Value Reduced Cost X1 10.00000 0.0000000 X2 10.00000 0.0000000 X3 8.000000 0.0000000 R11 3.000000 0.0000000 R12 2.000000 0.0000000 R13 0.0000000 0.0000000 R21 0.0000000 0.0000000 R22 1.000000 0.0000000 R23 0.0000000 0.0000000 R31 1.000000 0.0000000 R32 1.000000 0.0000000 R33 0.0000000 0.0000000 R41 0.0000000 0.0000000 R42 0.0000000 0.0000000 R43 2.000000 0.0000000

Row Slack or Surplus Dual Price 1 28.00000 0.0000000 2 0.0000000 0.0000000 3 0.0000000 0.0000000 4 0.0000000 0.0000000 5 1.000000 0.0000000 6 1.000000 0.0000000 7 0.0000000 0.0000000 8 3.000000 0.0000000 9 2.000000 0.0000000 10 3.000000 0.0000000 11 0.0000000 0.0000000 12 2.000000 0.0000000 13 3.000000 0.0000000

14 0.0000000 0.0000000 15 2.000000 0.0000000 程序8：

Max 0.1y1-0.2226x1+0.1833x2+0.2618x3+0.1695x4+0.1571y2+0.0196y3 st

1.5x1+2x2+x3+3x4<=14.4 x1+x2+x3<=5 x4<=2

y2-x1-2x2-4x4+y1=0

y3-10x1-4x2-16x3-5x4+2y1=0 y1-x1-2x2-4x4<=0

y1-5x1+2x2+8x3+2.5x4<=0 end

ROW SLACK OR SURPLUS DUAL PRICES

2) 0.000000 0.241577 3) 0.284615 0.000000 4) 0.000000 0.055237 5) 0.000000 0.157100 6) 0.000000 0.019600 7) 15.369231 0.000000 8) 0.000000 0.046373

NO. ITERATIONS= 3

RANGES IN WHICH THE BASIS IS UNCHANGED:

OBJ COEFFICIENT RANGES

VARIABLE CURRENT ALLOWABLE ALLOWABLE COEF INCREASE DECREASE Y1 -0.100000 0.342673 INFINITY X1 -0.222600 0.034507 1.570250 X2 0.183300 0.038297 0.028418 X3 0.261800 0.037162 INFINITY X4 0.169500 INFINITY 0.055237 Y2 0.157100 INFINITY 0.024155 Y3 0.019600 0.001725 0.078512

RIGHTHAND SIDE RANGES

ROW CURRENT ALLOWABLE ALLOWABLE RHS INCREASE DECREASE 2 14.400000 0.528571 6.900000 3 5.000000 INFINITY 0.284615 4 2.000000 1.840000 0.187342

5 0.000000 INFINITY 15.369231 6 0.000000 INFINITY 41.230770 7 0.000000 INFINITY 15.369231 8 0.00000

程序9：max 4.8x11+4.8x21+5.6x12+5.6x22-10x1-8x2-6x3 st

x-x1-x2-x3=0 x11+x12-x<500 x21+x22<1000 0.5x11-0.5x21>0 0.4x12-0.6x22>0 x1-500y1<0 x2-500y2<0 x3-500y3<0 x1-500y2>0 x2-500y3>0 end int y1 int y2 int y3

ENUMERATION COMPLETE. BRANCHES= 2 PIVOTS=

LAST INTEGER SOLUTION IS THE BEST FOUND RE-INSTALLING BEST SOLUTION...

OBJECTIVE FUNCTION VALUE

1) 5000.000

VARIABLE VALUE REDUCED COST Y1 1.000000 0.000000 Y2 1.000000 2000.000000 Y3 1.000000 1000.000000 X11 0.000000 1.200000 X21 0.000000 0.200000 X12 1500.000000 0.000000 X22 1000.000000 0.000000 X1 500.000000 0.000000 X2 500.000000 0.000000 X3 0.000000 0.000000 X 1000.000000 0.000000

ROW SLACK OR SURPLUS DUAL PRICES

28 2) 0.000000 6.000000 3) 0.000000 6.000000 4) 0.000000 5.000000 5) 0.000000 0.000000 6) 0.000000 -1.000000 7) 0.000000 0.000000 8) 0.000000 0.000000 9) 500.000000 0.000000 10) 0.000000 -4.000000 11) 0.000000 -2.000000

NO. ITERATIONS= 35

BRANCHES= 2 DETERM.= 1.000E 0

程序10：

Min 3x1+x2+3x3+3x4+x5+x6+3x7 s.t.

4x1+3x2+2x3+x4+x5>=50 x2+2x4+x5+3x6>=20 x3+x5+2x7>=15 end gin7

LP OPTIMUM FOUND AT STEP 4

OBJECTIVE FUNCTION VALUE

1) 26.66667

VARIABLE VALUE REDUCED COST X1 0.000000 1.666667 X2 11.666667 0.000000 X3 0.000000 1.666667 X4 0.000000 2.666667 X5 15.000000 0.000000 X6 0.000000 1.000000 X7 0.000000 1.666667

ROW SLACK OR SURPLUS DUAL PRICES 2) 0.000000 -0.333333 3) 6.666667 0.000000 4) 0.000000 -0.666667

NO. ITERATIONS= 4

RANGES IN WHICH THE BASIS IS UNCHANGED:

OBJ COEFFICIENT RANGES

VARIABLE CURRENT ALLOWABLE ALLOWABLE COEF INCREASE DECREASE X1 3.000000 INFINITY 1.666667 X2 1.000000 1.250000 1.000000 X3 3.000000 INFINITY 1.666667 X4 3.000000 INFINITY 2.666667 X5 1.000000 0.833333 0.666667 X6 1.000000 INFINITY 1.000000 X7 3.000000 INFINITY 1.666667

RIGHTHAND SIDE RANGES

ROW CURRENT ALLOWABLE ALLOWABLE RHS INCREASE DECREASE 2 50.000000 INFINITY 19.999998 3 20.000000 6.666667 INFINITY 4 15.000000 35.000000 9.999999

程序11： Model:

min=x1+x2+x3;

x1*r11+x2*r12+x3*r13>=50; x1*r21+x2*r22+x3*r23>=10; x1*r31+x2*r32+x3*r33>=20; x1*r41+x2*r42+x3*r43>=15; 4*r11+5*r21+6*r31+8*r41<=19; 4*r12+5*r22+6*r32+8*r42<=19; 4*r13+5*r23+6*r33+8*r43<=19; 4*r11+5*r21+6*r31+8*r41>=16; 4*r12+5*r22+6*r32+8*r42>=16; 4*r13+5*r23+6*r33+8*r43>=16; x1+x2+x3>=26; x1+x2+x3<=31; x1>=x2; x2>=x3;

@gin(x1);@gin(x2);@gin(x3); @gin(r11);@gin(r12);@gin(r13); @gin(r21);@gin(r22);@gin(r23); @gin(r31);@gin(r32);@gin(r33); @gin(r41);@gin(r42);@gin(r43); end

Rows= 15 Vars= 15 No. integer vars= 15

Nonlinear rows= 4 Nonlinear vars= 15 Nonlinear constraints= 4

Nonzeros= 73 Constraint nonz= 58 Density=0.304

No. < : 4 No. =: 0 No. > : 10, Obj=MIN Single cols= 0

Optimal solution found at step: 28970 Objective value: 28.00000 Branch count: 2946

Variable Value Reduced Cost X1 10.00000 0.0000000 X2 10.00000 0.0000000 X3 8.000000 0.0000000 R11 3.000000 0.0000000 R12 2.000000 0.0000000 R13 0.0000000 0.0000000 R21 0.0000000 0.0000000 R22 1.000000 0.0000000 R23 0.0000000 0.0000000 R31 1.000000 0.0000000 R32 1.000000 0.0000000 R33 0.0000000 0.0000000 R41 0.0000000 0.0000000 R42 0.0000000 0.0000000 R43 2.000000 0.0000000

Row Slack or Surplus Dual Price 1 28.00000 0.0000000 2 0.0000000 0.0000000 3 0.0000000 0.0000000 4 0.0000000 0.0000000 5 1.000000 0.0000000 6 1.000000 0.0000000 7 0.0000000 0.0000000 8 3.000000 0.0000000 9 2.000000 0.0000000 10 3.000000 0.0000000 11 0.0000000 0.0000000 12 2.000000 0.0000000 13 3.000000 0.0000000 14 0.0000000 0.0000000 15 2.000000 0.0000000 程序12：

Max 0.1y1-0.2226x1+0.1833x2+0.2618x3+0.1695x4+0.1571y2+0.0196y3 st

1.5x1+2x2+x3+3x4<=14.4 x1+x2+x3<=5 x4<=2

y2-x1-2x2-4x4+y1=0

y3-10x1-4x2-16x3-5x4+2y1=0 y1-x1-2x2-4x4<=0

y1-5x1+2x2+8x3+2.5x4<=0 end

ROW SLACK OR SURPLUS DUAL PRICES

2) 0.000000 0.241577 3) 0.284615 0.000000 4) 0.000000 0.055237 5) 0.000000 0.157100 6) 0.000000 0.019600 7) 15.369231 0.000000 8) 0.000000 0.046373

NO. ITERATIONS= 3

RANGES IN WHICH THE BASIS IS UNCHANGED:

OBJ COEFFICIENT RANGES

VARIABLE CURRENT ALLOWABLE ALLOWABLE COEF INCREASE DECREASE Y1 -0.100000 0.342673 INFINITY X1 -0.222600 0.034507 1.570250 X2 0.183300 0.038297 0.028418 X3 0.261800 0.037162 INFINITY X4 0.169500 INFINITY 0.055237 Y2 0.157100 INFINITY 0.024155 Y3 0.019600 0.001725 0.078512

RIGHTHAND SIDE RANGES

ROW CURRENT ALLOWABLE ALLOWABLE RHS INCREASE DECREASE 2 14.400000 0.528571 6.900000 3 5.000000 INFINITY 0.284615 4 2.000000 1.840000 0.187342 5 0.000000 INFINITY 15.369231 6 0.000000 INFINITY 41.230770 7 0.000 13.400000 7.400000

第八次 微分方程数值解

P120实验内容四：

1.a :

欧拉和R-K和精确解的比较 function ydot=el1(x,y) ydot=y+2*x;

x=0:0.1:1; n=length(x); y(1)=1;

for i=2:n

y(i)=y(i-1)+0.1*(y(i-1)+2*x(i-1)) end

y1=3*exp(x)-2*x-2;

[x,y2]=ode45('el1',[0:0.1:1],1)

plot(x,y,'g',x,y1,'r',x,y2,'b+')

dsolve('Dy=y+2*x','y(0)=1','x')

dsolve('x^2*D2y+x*Dy+(x^2-1/4)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x')

2: 1）C小题

function ydot=el3(x,y)

ydot=[y(2);-y(2)/x-(1-1/(4*x*x))*y(1)];

[x,y]=ode45('el3',[pi/2:pi/100:3*pi],[2,-2/pi]); x1=pi/2:pi/10:3*pi

y1=sin(x1).*(2*pi./x1).^(1/2); plot(x1,y1,'r*',x,y(:,1))

4. 单摆：

一根长为l的 (无弹性)细线上端固定，下端系着一个质量为m的小球。在重力的

作用下小球处于平衡位置。使小球偏离平衡位置一个小角度，然后让它自由，小球

就会沿着圆弧摆动。在不考虑空气阻力的情况下，小球将作周期一定的简谐运动.

(1)试建立模型来描述单摆运动规律(初始偏角分别为

5??0.0873,30??0.5236).

(2)用matlab编写求解程序.(取l?25cm,g?9.8m?s?2) (1)分析与建立模型：

如图，以??0为平衡位置，以右边为正方向建立摆角?的坐标系.在小球摆动过程中的任一位置?,小球所受重力沿运动轨迹方向的分力为?mgsin?(负号表示力的方向与?的正向相反)，

本文来源：https://www.bwwdw.com/article/x0z6.html