最新数学建模第三次作业

院系：数学学院

专业：信息与计算科学年级：2014级

学生姓名：王继禹

学号：201401050335 教师姓名：徐霞

6.6 习题

3.一个慢跑者在平面上沿着他喜欢的路径跑步，突然一只狗攻击他，这只狗以恒定速率跑向慢跑者，狗的运动方向始终指向慢跑者，计算并画出狗的轨迹。

解：

（1）模型分析建立:

狗的轨迹：在任意时刻，狗的速度向量都指向它的目标慢跑者。

假设1：慢跑者在某路径上跑步，他的运动由两个函数X(t)和Y(t)描述。

假设2：当t=0时，狗是在点(x0,y0)处，在时刻t时，它的位置是(x(t),y(t))

那么下列方程成立：

(1)狗以恒定速率跑: X’2+y’2=w2

(2) 狗的速度向量平行于慢跑者与狗的位置的差向量：

将上述方程带入等式：,可得：

再将λ代入第二个方程，可得狗的轨迹的微分方程：

（2）程序及结果

dog函数

[dog.m]

function [zs,isterminal,direction] = dog(t,z,flag)

global w;% w=speed of the dog

X=jogger(t);

h = X-z;

nh=norm(h);

if nargin<3 || isempty(flag)

zs=(w/nh)*h;

else

switch(flag)

case'events'

zs = nh-1e-3;

isterminal = 1;

direction = 0;

otherwise

error(['Unknow flag:' flag]);

end

慢跑者的运动轨迹方程，水平向右

[jogger.m]

function s = jogger(t);

s = [8*t;0];

标记的函数

[cross.m]

function cross(Cx,Cy,v)

Kx = [Cx Cx Cx Cx-v Cx+v];

Ky = [Cy Cy+2.5*v Cy+1.5*v Cy+1.5*v Cy+1.5*v] plot(Kx,Ky);

plot(Cx,Cy,'o');

主程序：

静态显示

[main1.m]

global w

y0 = [60;70];

w=10;

options = odeset('RelTol',1e-5,'Events','on'); [t,Y] = ode23('dog',[0,20],y0,options);

clf;hold on;

axis([-10,100,-10,70]);

plot(Y(:,1),Y(:,2));

J=[];

for h=1:length(t),

w = jogger(t(h));

J=[J;w'];

end

plot(J(:,1),J(:,2),':');

p = max(size(Y));

cross(Y(p,1),Y(p,2),2)

hold off;

动态显示

[main2.m]

global w;

y0=[60;70];

w=10;

options = odeset('RelTol',1e-5,'Events','on'); [t,Y]=ode23('dog',[0,20],y0,options); J=[]; for h=1:length(t);

w= jogger(t(h));

J=[J;w'];

xmin = min(min(Y(:,1)),min(J(:,1)));

xmax = max(max(Y(:,1)),max(J(:,1)));

ymin = min(min(Y(:,2)),min(J(:,2)));

ymax = max(max(Y(:,2)),max(J(:,2)));

clf;hold on;

axis([xmin-10 xmax ymin-10 ymax]);

title('The jogger and the Dog');

for h = 1:length(t)-1,

plot([Y(h,1),Y(h+1,1)],[Y(h,2),Y(h+1,2)],'-','Color','red','EraseMode ','none');

plot([J(h,1),J(h+1,1)],[J(h,2),J(h+1,2)],'-','Color','green','EraseMo de','none');

drawnow;

pause(0.1);

end

plot(J(:,1),J(:,2),':');

p = max(size(Y));

cross(Y(p,1),Y(p,2),2)

hold off;

结果

t=12.2761812635281，在12.27秒后狗追上慢跑者。

慢跑者轨迹是椭圆轨迹

[jogger2.m]

function s=jogger2(t)

s=[10+20*cos(t) 20+15*sin(t)];

狗的微分方程

[dog.m]

function [zs,isterminal,direction] = dog(t,z,flag) global w;% w=speed of the dog

X=jogger2(t);

h = X-z;

nh=norm(h);

if nargin<3 || isempty(flag)

zs=(w/nh)*h;

else

switch(flag)

case'events'

zs = nh-1e-3;

isterminal = 1;

direction = 0;

otherwise

error(['Unknow flag:' flag]);

end

end

主程序

[main3.m]

global w;

y0=[60;70];

w=10;

options = odeset('RelTol',1e-5,'Events','on');

[t,Y]=ode23('dog',[0,20],y0,options); J=[];

for h=1:length(t);

w= jogger2(t(h));

J=[J;w'];

end

xmin = min(min(Y(:,1)),min(J(:,1)));

xmax = max(max(Y(:,1)),max(J(:,1)));

ymin = min(min(Y(:,2)),min(J(:,2)));

ymax = max(max(Y(:,2)),max(J(:,2)));

clf;hold on;

axis([xmin-10 xmax ymin-10 ymax]);

title('The jogger and the Dog');

for h = 1:length(t)-1,

plot([Y(h,1),Y(h+1,1)],[Y(h,2),Y(h+1,2)],'-','Color','red','EraseMode ','none');

plot([J(h,1),J(h+1,1)],[J(h,2),J(h+1,2)],'-','Color','green','EraseMo de','none');

drawnow;

pause(0.1);

end

plot(J(:,1),J(:,2),':');

p = max(size(Y));

cross(Y(p,1),Y(p,2),2)

hold off;

结果取w=25有

t=4.017776368842910，经过4秒左右狗追上慢跑者。

8.平面上有n(n>=2)个圆，任何两个圆都相交但无3个圆共点。试问n个圆把平面划分成多少个不连通的区域？

解：∵一个圆将平面分为2份

两个圆相交将平面分为4=2+2份，

三个圆相交将平面分为8=2+2+4份，

四个圆相交将平面分为14=2+2+4+6份，

…

平面内n个圆，其中每两个圆都相交于两点，且任意三个圆不相交于同一点，

则该n个圆分平面区域数f（n）=2+（n-1）n=n2-n+2

证明：（1）当n=1时，一个圆把平面分成两个区域，而12-1+2=2，命题成立．

（2）假设n=k（k≥1）时，命题成立，即k个圆把平面分成k2-k+2个区域．

当n=k+1时，第k+1个圆与原有的k个圆有2k个交点，这些交点把第k+1个圆分成了2k段弧，

而其中的每一段弧都把它所在的区域分成了两部分，因此增加了2k个区域，

共有k2-k+2+2k=（k+1）2-（k+1）+2个区域．

∴n=k+1时，命题也成立．

由（1）、（2）知，对任意的n∈N*，命题都成立．

9.某人有n元钱，他每天买一次物品，每次买物品的品种很单调，或者买一元钱的甲物品，或者买二元钱的乙物品，问他花完这n元钱有多少不同的方式？

解：设an表示花完这n元钱的方案种数，

若n=1，则只能买甲，有一种方法，故a1=1，

若n=2，则可以买2个甲，或者1个乙或1个丙，即a2=3，

当n≥3时，花钱的方式由购买甲和购买乙购买丙的种数之和构成，

即a n=a n-1+a n-2+a n-2=a n-1+2a n-2

则当n≥3时，a n+a n-1=2（a n-1+a n-2），

即{a n+1+a n}是公比q=2的等比数列，首项为a2+a1=1+3=4，

则a n+1+a n=4?2n-1=2n+1，

∴a n+a n-1=2n，

两式相减得a n+1-a n-1=2n+1-2-=2-，（n≥2），

若n是奇数，a n=2n-1+2n-3+…+22+a1=(2n+1-1)/3

若n是偶数，a n=2n-1+2n-3+…+23+a2=(2n+1+1)/3．

7.6 习题

1.在化工生产中常常需要知道丙烷在各种温度T和压力P下的导热系数K。下面是实验得到的一组数据：

10KPa下的K。

试求T=99C?和P=10.3x3

解：找出温度T相等时，导热系数K与压力P的关系。由于在不同温度时，仅给出两个K、P的值，因此采用线性近似，把K、P看作是线性关系。

建立M文件：

function y=y_lagr1(x0,y0,x)

n=length(x0);m=length(x);

for i=1:m

z=x(i);

s=0.0;

for k=1:n

p=1.0;

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

主程序：

p1=[9.7981,13.324]; k1=[0.0848,0.0897]; %T=68℃

p2=[9.0078,13.355]; k2=[0.0762,0.0807]; %T=87℃

p3=[9.7918,14.277]; k3=[0.0696,0.0753]; %T=106℃

p4=[9.6563,12.463]; k4=[0.0611,0.0651]; %T=140℃

a2=polyfit(p2,k2,1); a3=polyfit(p3,k3,1);

x1=polyval(a2,10.3); x2=polyval(a3,10.3);

%x1，x2分别是P=10.3*10^3kPa下87℃和106℃的K值

plot(10.3,x1,'k+',10.3,x2,'k+',p1,k1,p2,k2,p3,k3,p4,k4)

xlabel('丙烷压力P')

ylabel('丙烷导热系数K')

title('在不同温度下丙烷导热系数与压力的关系图')

gtext('T=68℃'),gtext('T=87℃'),gtext('T=106℃'),gtext('T=140℃')运行后图中所标点为P=10.3*10^3kPa时，T=87℃和T=106℃对应的导热系数K值。

在T=87℃和T=106℃之间仍采用线性近似来求T=99℃时的导热系数K。

程序如下：

x=[87,106];

y=[x1,x2];

a=polyfit(x,y,1);

z=polyval(a,99)

z=0.0729

plot(99,z,'k+',x,y)

grid

xlabel('丙烷温度T')

ylabel('丙烷导热系数K')

title('压力P=10.3*10^3kPa时丙烷导热系数与温度的关系') 运行结果: T=99℃、P=10.3*10^3KPa时K=0.0729。

4.用电压V=10伏的电池给电容器充电，电容器上t时刻的电压为v(t)=V-(V-V0)

t

e T-,其中v0

解：

建立M文件function f=dianya(x,t)

f=10-(10-x(1))*exp(-t/x(2))

%x(1)=V0；x(2)=τ

主程序：

t=[0.5 1 2 3 4 5 7 9];

V=[6.36 6.48 7.26 8.22 8.66 8.99 9.43 9.63];

x0=[0.2,0.05];

x=lsqcurvefit('dianya',x0,t,V) %x(1)代表V0，x(2)代表τ

f=dianya(x,t)

结果

x=5.5577 3.5002,即初始电压v0=5.5577，充电常数为τ=3.5002．