2023年MATLAB数学实验报告

MATLAB数学实验报告

实验日期：2023.5

学院：能源与动力工程

班级：化工11

成员：王旭陆清华仲秋晨

一、实验目的

1.学习MATLAB软件的循环和选择结构，进一步提高MATLAB编程能力；

2.通过对一些基础数学实验的学习和实践，了解级数逼近和数值积分、用

最小二乘法进行数据拟合等的数学思想和数学方法，开拓数学视野，提

高数学水平。

二、实验内容

1.⑴问题:对于数列{Jn},n=1,2,・・•，求当其前n项和不超过1000

时的值以及和的大小。

(2)分析：这个问题书上已有例题解答，但是书上的程序运营结果

最后一行结果并不是我们所规定的解的答案，而倒数第二行则是所求问

题的解。以下是修改后的程序已解决此问题。

(3)程序：

cIear;clc;

n=1;

s=1;

whilesV=1000

fprintf('n=%.Of,s=%.4f\n',n,s)

n=n+1;

s=s+sqrt(n);

end

(4)运营结果

n=123,s=914.7651

n=124,s=925.9007

n=125,s=937.0810

n=126,s=948.3060

n=127,s=959.5754

n=128,s=970.8891

n=129,s=982.2469

n=130,s=993.6487

2.(1)问题：1790年到1980年各年美国人口数的记录数据如下表:

美国人口记录数字(单位：百万)

年

1718001810182018301818186018701880

份

904050

记3.95.37.29.612.917.123.231.438.650.2

录

年

189011910192019301941950191971980

份

9000600

记62.72.092.106.5123.213150.717204.022

录001.79.36.5

是根据前12023的数据，分别用MaIthas模型和Logistic模型速立美国

人口增长的近似曲线(设美国人口总容纳量为1。亿)，并预测后2023的人口数，

通过与实际数据相比较,对两种预测结果进行分析。

(2)分析:根据题目规定分别用MaIthas模型和Logistic模型建立美国人口

增长的近似曲线。

(3)程序：

%MaIthas

cIear;clf

t=1790：10：1980；

N=[3.95.37.29.612.917.123.231.438.650.262.072.0

92.0106.5123.2131.7150.7179.3204.0226.5];

pIot(t,N,k.','markersize',20);

axis([179020803400])；

grid；hoIdon

pause(0.5)

n=20；

a=sum(t(1：n))；

b=sum(t(1：n).*t(1：n))；

c=sum(log(N(1:n)));

d=sum(t(1:n).*Iog(N(1：n)))；

A=[na;ab]；

B=[c；d];

p=inv(A)*B

x=179O:2023;

y=exp(p(1)+p(2)*x)；

pIot(x,y,'r—','Iinewidth',2)

%Logistic

cIear;cIf

t=179O:10:1980;

N=[3.95.37.29.612.917.123.231.438.650.262.07

2.092.0106.5123.2131.7150.7179.3204.0226.5];

pIot(t,N,'k.','markersize',20);

axis([1790208001000]);

grid：hoIdon

pause(0.5)

n=20;

k=1000;

N<-1-kf

a=sum(t(1:n));

b=sum(t(1:n).*t(1:n));

c=sum(Iog(M(1:n)));

d=sum(t(1:n).*log(M(1：n)))；

A=[na;ab];

B=[c;d];

p二inv(A)*B

x=1790：10：2080；

y=/k)+exp(p(1)+p(2)*x));

pIot(x,y,'r-',Iinewidth',2)

3.(1)追击问题:在一边长为1的正方形跑到的四个顶点上各站有1人,

他们同时开始以等速度沿跑道追逐下一个人,在追击过程中，每个人时

刻对准目的，试模拟追击路线。

(2)分析：参考书上导弹追击飞机问题。设4人为A,B,C,D,则先把A看做

飞机，把追击A的B看做导弹进行模拟；再把B看做飞机，把追击B的C看做导

弹进行模拟，以此类推。

(3)程序：

clear;cIc;cIf;

holdon

axis([O1100110]);

grid

A=[0,0];

B=[100,0];

C=[100,100]；

D=[0,100]；

k=0；

v=1；dt=1;

whiIek<10000；

k=k+1；

pIot(A(1),A(2),'r.','markersize',15);

plot(B(1),B(2),'b.','markersize',15);

pIot(C(1),C(2),y.','markersize',15);

plot(D(1),D(2),'g.','markersize',15)；

d=norm(A-B);e=(A-B)/d;

fprintf('k=%.0fA(%.2f,2f)B(%.2f,%.2f)C(%.2「％.2f)

D(%.2f,%.2f)d=%.2f\n',k,A(1),A(2),B(1),B(2),C(1),0(2),D(1),

D(2),d)

ifd<1.5

break

end

B=B+v*dt*e；

d=norm(B-C);e=(B-C)/d;

C=C+v*dt*e;

d=norm(C-D);e=(C-D)/d;

D=D+v*dt*e;

d二norm(D—A)；e=(D—A)/d;

A=A+v*dt*e；

pause(0.2)

end

(4)结果

4CommandWindow

fileEditDefaugDesktopJtfndowHelp

1•'NewtoMATLAB7WatdithisVideo,seeDemos,orreadGettingStarted.

k=70AGO.05,5S.02)B(41.64,29.43)CCO.05,41,87)D(58.03,69.96)4=30.55

k=71A(30.98,58.39lB(41.26,30.36)C(69.12,4!.5O)D(58.40,69.03)d=29.85

k=72A⑶.92,58.：2)8(40,92,31.30)C(68.18,41.16)M限73,68.09)*28.86

k=73A(32.88,59.0115(40.61,32,25)C(67.23,40.$5)D(59.02,67.13)4:2：.86

k=74A(33.84,59.27)5(40.33,33.21)C(66.2*,40.58)D(59.23,66.17)d:26.86

k=75A04.82,59.50)8(40.09,34.18)C(65.30,40.34)D(59.5】,65.19)*25.86

k=：6A(35.80,59.68)6(39.88,35.16)€(64.32,40.14)D(59.70,64.21)4:24.$6

k=77A(36.79,59.83)8(39.72,36.15)063.33,39.98)D(59.85,63.22)d=23.8；

k=78A(37.78,69.93)8(39.60,37.14)C(62.34,39.86)D(59.96,62.23)4:22,87

k=79A(38.78,59.9916(39.52,38.13)C(61.34,39.1)D(60.8,61.23)£21.8；

k=80A(39.78,60.00l5(39.48,39.13)C(60.34,39,75)D(60.03,6。.23)£20.8：

k=81A(40.78,59.98)6(39.50,40.C(59.34,39.77)D(M.00,59.23)4:19,87

k=82A«1.78,59.87)8(39.56,41.13)0(58.34,39.84)D(59.91,58.23)d=18,87

k=S3A(42.7；59.：3)B(39.68,42/2)C(5；.35,39.96)D(59.77,57.24)4=17,87

k=84A(43.75,59.53)5(39.85,43.11)C(56.37,40.14)D(59.58,56.26)*16.8；

k=85A(«.Tl,59.27)8(40.0834.08)0(55.39,W.38)D(59.瑞55.29)*15.8：

k=86A(45.66,58.94}B(40,37,45.04)C(54.U,40.67)M59.01,54.35)*487

k=87A(46.58,53.55)B(40.73,45.97)C(53.51,41.03)D(58.63,53.42)d=13,8T

k88A街.隔58.09)B(41.15,46.88)C(52.60,41.46)D(58.18,52.53)4=12,86

k=89A(48.3l,57,56)BUl.64J7.75)C(5d.96)D(57.66,51.67)d=ll.86

k=90A(49.10,56.95l5(42.20,48.58)0(50.91,42,53)D(57.06,50.87)4:10.$5

k=91A(49.83,M.26)班42.84,49.35)C(50.15,4117)D(56.39,50.13)29.83

lt=92A(50.48,55.50)B«3.65,S0.05)C(49.46J3.90)0(55.68,49.46)^8.82

k=93A(51.03,54.67)5(44.54,50.67)C(48.85J4,69)D(54.83,48.88)d=7.80

k=94A(51.46,53.76)5(45.20,51.18)0(48.36,45.57)D(53.94,48.43)d=6.77

k=95A(51.72,52.80)8(46.12,51.57)CG8.0lJ6.50)0(62.99,48.12)d：5.73

k=96A(51.T8,51.80l6(47.10,51.78)CU7.8437.49)9(52.00,48.00)^4.680申^^”匕|

k=9：A(51.57,50.82)B(48.10,51.78)€(47.90,48.49)D(51.00,48.12)4:3.60

k=98A(51.03J9.98)8(49.06,51.52)C(48.26,49.42)D(50.10,48.50*2.50

k=99A(50.10J9.60)8(49.85,50.90)CM8.99,50,10)D(帆52,49.36)*1.33

k»

舌♦耶.际e

4.(1)问题：【水塔流量的估计】美国某州的用水管理机构规定各社区提

供以每小时多少加仑计的用水量以及天天的用水量。许多社区没有测量

流入或流出水塔水量的装置，只能代之以每小时测量水塔中的水位，其

误差不超过5%。淡水他天天有1到2次的水泵供水，每次约两小时。当水

塔中的水位下降到最低水位L时，水泵就自动向水塔输水直到最高水位

H,此期间不能测量水位。现在，已知该塔是一个高40ft(英尺)，直径

57代的正圆柱，某小镇一天水塔水位的记录数据为下表：

某小镇某天水塔水位记录

时间(S)水位(ft)时间(s)水位(ft)

031.754663633.50

331631.104995332.60

663530.545393631.67

1061929.945725430.87

1393729.476057430.12

1792128.926455429.27

2124028.506853528.42

2522327.957185427.67

2854327.527502126.97

3228426.9779254水泵启动

35932水泵启动82649水泵启动

39332水泵启动8596834.75

3943535.508995333.97

4331834.459327033.40

其中水位降至约27ft水泵开始工作，水位上升到35.5ft时停止工作。

试估计任何时刻t（涉及书泵工作时间）从水塔中流出的水流量Q（t）,并

估计一天的总水量。

（2）分析：

1）拟合水位-时间函数

测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段），

和3个水泵不工作时段（以下称第1时段t=0到t=32284,第2次时段t=39435

到t=75O21和第3时段t=85968以后）.对第1、2时段的测量数据直接分

别作多项式拟合，得到水位函数.为使拟合曲线比较光滑，多项式次数不要太高,

一般在3〜6.由于第3时段只有3个测量记录，无法对这一时段的水位作出较好

的拟合.

2）拟定流量~时间函数

对于第1、2时段只需将水位函数求导数即可，对于两个供水时段的流量,

则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2

供水时段流量外推，将第3时段流量包含在第2供水时段内.

3）一天总用水量的估计

总用水量等于两个水泵不工作时段和两个供水时段用水量

之和，它们都可以由流量对时间的积分得到。

（3）程序：

t1=[033166635106191393717921212402522328543322

843943543318466364995353936572546057464554685

35718547502185968899539327O];

h=[31.7531.1030.5429.9429.4728.9228.5027.9527.52

26.9735.5034.5033.5032.6031.6730.8730.1229.2728.42

27.6726.9734.7533.9733.40]；

t=t1./360O;

s=pi*(57*0.3048/2)-2；

v=h.*O.3048.*s.*1000./3.78541

c1=poIyfit(t(1:10),v(1:10),3);

a1=poIyder(d);

tp1=0：0.1:9;

x1=-polyvaI(a1,tp1)；

fprintf('f1(t)=%.4fx"2+%.4fx+%.4f\n',a1(1),a1(2),a1⑶)；

plot(tp1,x1,,g-')

c2=polyfit(t(11：21),v(11：21),3)；

a2=poIyder(c2);

tp2=11：0.1:21;

x2=-poIyvaI(a2,tp2);

fprintfCf3(t)=%.4fx^2+%.4fx+%.4f\na2(1),a2(2),a2(3))；

hoIdon

pIot(tp2,x2,'b-1)

xx1=-poIyvaI(a1,[89])；

xx2=-polyval(a2,[1112]);

xx12=[xx1xx2]；

d2=polyfit([891112],xx12,3);

tp12=9：0.1：11；

x12=poIyval(c12,tp12);

fprintf('f2(t)=%.4fxe3+%.4fxe2+%.4fx+%.4f\n',c12(1),c

12(2),c12(3),c12(4))；

holdon

pIot(tp12,x12,'y-')

dt3=diff(t(22：24))；

dv3=diff(v(22：24))；

dvt3--dv3./dt3;

t3=[2021t(22)t(23)];

a-[~Polyval(a2,t3(1：2))]；

xx3=[advt3];

c3=polyfit(t3,xx3,3)；

tp3=21:0.1:24;

x3=poIyvaI(c3,tp3);

fprintf('f4(t)=%.4fx八3+%.4fx人2+%.4fx+%.4f\n',c3(1),

c3⑵，c3(3),c3(4));

holdon

plot(tp3,x3,'r-')

y1=0.1trapz(x1);

y2=0.1*trapz(x2)；

y12=0.1*trapz(x12);

y3=0.1*trapz(x3)；

y=(y1+y2+y12+y3)

grid

(4)结果：

f1(t)=-143.7890x^2+1656.5813x+-13731.5134

f3(t)=-86.8481x^2+3341.3267x+-46996.7559

f2(t)=-903.1315xA3+26738.6695x^2+-257790.0151x

+823129.9568

f4(t)=-3.8291x^3+41.6810x^2+3339.0536x+-37911.0

031

总水量y=

3.3167e+05

5.(1)问题:试按照级数逼近和数值积分的思想方法，计算n的近似值(精

确到10、5)。规定：

1)运用级数展开公式(7-2)—(7-6)来计算。

2)运用梯形法公式(7-1。)、抛物线法公式(7T2)分别进行计

算并加以比较。

(2)程序:

1)7-2

cic;cIear；

n=0;r=1;

p=0;k=-1;

whiIer>=1.Oe-5

n=n+1;

k=k*(—1);

p1-p+k/(2*n-1);

r=abs(4*(p1-p));

fprintf('n=%.Of,p=%.10f\n',n,4*p1)；

P=p1；

end

运营结果：

n=199979,p=3.

n=199980,p=3.

n=199981,p=3.

n=199982,p=3.

n=199983,p=3.

n=199984,p=3.

n=199985,p=3.

n=199986,p=3.

n=199987,p=3.

n=199988,p=3.

n=199989,p=3.

n=199990,p=3.

n=199991,p=3.

n=199992,p=3.

n=199993,p=3.

n=199994,p=3.

n=199995,p=3.

n=199996,p=3.

n=199997,p=3.

n=199998,p=3.

n=199999,p=3.

n=202300,p—3.

n=202301,p=3.

7-3

cic；clear；

n=0；r=1;p=0;

k=-1;a=1；b=1;

whiler>=1.Oe—5

n=n+1;

k-k*(-1);

a-4*a；b=9*b；

p1=p+k/(2*n-1)*(2/a+3/b);

r=abs(4*(p1-p))；

fprintfCn=%.Of,p=%.10f\n',n,4*p1)；

P=p1；

end

运营结果：

n—1,p—3.

n—2,p—3.

n=3,p=3・

n=4,p=3.

n--5,p—3.

n—6,p=3.

n--7,p=3.

n=8,p=3.

7-4

clc;clear；

n=0;r=1;p=0;

k=-1；a=1;b=1;

whiler>=1.Oe—5

n=n+1;

k=k*(-1)；

a=25*a；b=239*239*b;

p1=p+k/(2*n-1)*(20/a-239/b)；

r=abs(4*(p1-p))；

fprintf(fn=%.0f,p=%.10f\n',n,4*p1)

p=p1;

end

运营结果：

n—1,p—3.

n=2fp=3.

门=3,p=3.

n=4,p=3.

n=5,p-3.

7-5

clc;cIear;

n=0；r=1；p=0；

k=-1；a=1;b=1；

whiIer>=1.0e-5

n=n+1;

k=k*(T);

a=9*a;b=49*b；

p1=p+k/(2*n-1)*(6/a+7/b)；

r=abs(4*(p1-p));

fprintf('n=%.Of,p=%.10f\n',n,4*p1)；

P=p1；

end

运营结果：

n=1,p=3.

n—2,p—3.

n—3,p=3.

n=4,p=3.

n=5,p=3.

n=6,p—3.

»

7-6

cic;clear;

n=0;r=1;p=0;k=-1;

a-1;b=1;c=1;d=1;

whiler>=1.0e-5

n=n+1;

k=k*(-1)；

a=28*28*a;b=443*443*b；c=1393*1393*c；d=11018*1101

8*d；

p1=p+k/(2*n-1)*(22*28/a+886/b-