数学模型与数学建模实验三

实验报告三

学院名称：理学院 专业年级： 姓 名： 学 号：

课 程：数学模型与数学建模 报告日期：2015年11月24日

一、实验题目

例1.3.1 继续研究十字路口红绿灯问题 十字路口绿灯亮灯30s ，最多可以通过多少辆汽车？继续研究问题：十字路口绿灯亮t s ，最多可以通过多少辆汽车)(t n n =？ 例 1.3.2 图像做旧

练习：改变融合比例，尝试其它做旧效果 例1.3.3 画分形树

二、实验目的

数学建模是一种数学的思考方法，用数学的语言和方法，通过抽象、简化建立能近似的刻画并“解决”实际问题的路径。因此，模型一般具有明确的应用背景，针对性较强。为了使所研究的模型有相当的普适性，能解决一类问题，就需要在模型确立之后，进一步分析推广，挖掘出模型更多的理论和实用的价值。

三、问题陈述

第一题

问题要求很明确，求解的关键是绿灯亮后，汽车才启动，加速是过路口，要确定在给定时间内能通过多少辆车，就要考虑汽车加速度，从停车位置到路口的路程以及城市行车的最高限速要求。 第二题

数值图像是一个函数，4:R f →Ω，确切的说，是定义在矩形区域3

R →Ω的离散网格点()i i y x ,上的函数，称()i i y x ,为像素点()i i ij y x f z ,=为像素值，当1=d 时，f 为灰度图像，当3=d 时，f 为彩色图像，所以灰度图像就是一个矩阵()

ij z Z =，彩色图像就是三个矩阵，分别表示三个颜色（红、绿、蓝RGB ），即三维数组表达，图像做旧是一种图像融合的方法，可以通过两个矩阵的加权求和来实现。

需要强调的是，为解决问题编写程序不仅需要熟悉软件的变异特点，例如MATLAB 的矩阵语言，而且要密切结合问题的背景，利用其数学特征，才有可能编写出有效的程序，下面以编写分形树的绘图为例，说明如何利用分形树的自相似特征，将需要运行十多分时间的程序改进为几秒就能实现的程序。

第三题

分形是对不规则的难以用传统欧式几何描述的几何图形，例如，海岸线和山川形状，多数分形图案的特点之一是从整体到局部的自相似性，从远距离观察，海岸线和山川形状是极不规则的；从近距离观察，其局部形状又和整体形态相似，数学的分形树是按一定规律构造

的的具有自相似特征的几何图案。我们从一个线段生成一个分形树图案开始，学习几何最基本的元素“点”和“直线”的矩阵代数表达方式。

四、模型及求解结果

第一题

记在时刻t 第n 辆车的位置为)(t S n ，用数轴表示车辆行驶道路，数轴的正向为汽车行驶方向，数轴原点为红绿灯的位置。于是当0)30(>n S 时，表明在第30s 第n 辆车已通过红绿灯，否则，结论相反。于是，只要确定限速行驶模型)(t S n ，就可以确定30s 内通过的汽车数量n 。在单向、单车道、直行、限速等假设下，在《数学模型与数学建模》中给出了这个模型：

???

??≤≤≤≤-+≥-+-+=.

0),0(,

,)(2/)0(,),()(2/)0()(*2**2*n n

n n n n n n n n n n t t S t t t t t a S t t t t v t t a S t S 其中))(1()0(D L n S n +--=表示第n 辆车在绿灯亮前的位置，n t n =表示第n 辆车的启动时间，n n t a v t +=/**表示第n 辆车达到最大限速时间。参数驱车长m L 5=，车距m D 2=，加速度2

/2s m a =，最大限速2

*/11s m v =。可以得到结果，绿灯亮30s ，该路口单向，单车道可以通过17辆车。

这里直接研究限速行驶模型的应用。确定十字路口绿灯亮多长时间是城市交通管理中最基本的一个问题。直接利用限速行驶模型计算每辆车到达路口的时间，有数值结果可见，绿灯亮后汽车开始以最高限速穿过路口的时间在11s 以后，从第6辆车开始。以后的车都以最高限速穿过路口。绿灯亮30s 通过路口17辆车中有11辆车以最高限速穿过路口，如果绿灯只亮15s ，则在通过7辆车中只有两辆车以最高的限速穿过路口，显然这样的交通灯控制策略对于路口的利用率是不高的。

如果是利用限速行驶模型得到通过的车辆数对绿灯亮的时间长度的依赖关系)(t n n =，即得到当绿灯亮t s 时，单向、单车道通过的车辆数，对进一步的研究更方便。要写出这个分段函数，推导并不困难，但是有些繁琐。但是在实际问题研究中找到显示表达式函数的机会非常少，借助计算机编程，；iyong 输入输出表达函数关系常常是唯一可取的办法。下面以这个简单问题为例尝试一下。

一个基本的想法是，在有限时间内只能通过有限辆车，只要计算出n 辆车通过路口的时间，再用规定的绿灯亮的时间长度为卡尺，就可以确定通过的车辆数。 算法

（1）计算每辆汽车达到最大限速时间0t ，以加速度通过路口时间1t 和以最大限速通过路口时间2t 。

（2）比较0t 和1t 就可以确定这辆车实际通过路口的时间。

（3）将充分多辆车通过路口的时间列出来，（用find ）确定通过路口时间小于绿灯亮的时间长度T 的所有车号，其中最大者就是能够通过的车辆数，函数以pass.m 文件名存储。 于是，只要输入T 就会由函数pass(T)得到通过的车辆数。 第二题

1、图像做旧处理所得图像及其拉伸灰度级显示图像结果如下图所示：

上图 （图像做旧）左上为原双精度图，右上为整数型数据的得到的图像，左下为双精度整数型数据各50%的融合比例，右下80%的双精度数据与20%的整数型数据融合后的结果。 2、练习题：改变融合比例

上图 （改变融合比例）左上为原双精度图，右上为整数型数据的得到的图像，左下为40%的双精度数据与60%的整数型数据融合后的结果，右下90%的双精度数据与10%的整数型数据融合后的结果。

第三题

第（1）步：画一个点P （2.5,3）

第（2）步：画两个点()3,5.21P

,()6,5.62P

第（3）步：画两个点1P

，2P 及其连线上的中点，因为任意两点()111,y x P ，()221,y x P 之间的连线上的点

()αααy x P ,可以表示为

()211P P P ααα-+=，即：()()??

?-+=-+=212

111y y y x x x αααααα

其中10≤≤α，特别地，当2/1=α时，αP 是连线上的中点。

注意：两个同阶矩阵相加等于对应元素相加，矩阵除以某数等于每个元素除以这 个数。

第（4）步：一枝树杈，任意给定一节树段，在其中点长出一个分杈，长度等于原树段长的一半，向左偏离原树段30°角，记原树段起点为1P ，终点为2P ，中点为0P ，如图：

算法：

①将向量21P P

平移到原点，并缩小1/2，得到向量OP

，

()2

/12P P P -=

②将

OP

逆时针旋转π/6，再平移到1P 与2P 连线的中点2102/12/1P P P +=，得到向量

L P P 0，PA P P L +=0,其中

?

??? ??-=6/cos 6/sin 6/sin 6/cos ππππA 是旋转变换矩阵；

③将1P

，0P ，L P ，0P ，2P 依序连接起来。

第（5）步：画一棵分形树。假定树的生长是有规律的，按分杈点可以将树分成有限个树段，每个阶段每个树段必在且只在中点向一个固定方向，例如，向左偏离30°，长出一枝新树杈，新树杈的出现，使得树可以被分为更多的树段，进入下一阶段的生长，在计算机上实现这个过程，上一步已将一个树段变为一个树杈，这树杈由3个新树段构成，为表达清楚，可将上题的的输出指令

u=[u(1,:);p0;Ip;p0;u(2,:)] 改为

u=[u(1,:);p0;p0;lp;p0;u(2,:)]

即将每个树段的起点和终点坐标都分别用矩阵u 的奇数行和偶数行记录，3个树段由6个点表示。下一步，在这3个树段的每段上都要长出一个分杈，即在每段上重复上一步工作，为编写代码简洁，先将画一枝树杈的指令用函数文件branch.m 表达。

只要输入树段断电坐标v 和将要长出分杈的偏离角度θ，就会输出长了一个分杈的树杈，然后，重复执行调用branch.m 画出所有树段生成的树杈，这个过程用函数文件tree.m 表达

只要输入上阶段已有的树u 和规定的树杈偏离角θ，就会输出长了更多树杈的树。注意到，因为u 记录了所有树段的起点和终点，所以树段总数是其行数的一半。按这种想法，可以设计生长0k 个阶段的分形树的算法。

算法

①:0=k 给定初始值0u ，一个2行2列矩阵表示一根树干两端点的坐标，给定参数θ表示

树杈的生长角度。 ②第k 阶段（1≤k ）： a.从上阶段的树

1-k u 截出树段l v ；

b.树段l

v 上长成一枝树杈，记为l

uu ；

c.依序将l

uu 拼接，构成k u

③如果k <

0k ,则k =k +1，重复②，否则结束运算。

以下分别画出k =0,1，...，8阶段的分形树。

如图，注意到，不论k 多大，也就是不管迭代多少次，树杈越来越密，但是永远不会充满整个区域，这个特点可以利用Matlab 同性工具栏的放大镜看到。于是，这棵树的树杈长度综合可无限增大，但是树杈占据的面积永远为零，这就是所谓的二维分形。分形的一个重要特征就是自相似，即局部是整体的相似形，通过图形放大镜可以展现这个特征。标准的自相似分形是数学对现实的抽象，实际分形的形体自相似性大多是统计意义上的。

如果画一颗每阶段在每段数段上长两个分叉的分形树，一个分叉向左，一个分叉向右，就会发现按照上面的思路够早的算法计算速度很慢。可以在程序开始前加上指令tic ，在结尾后加上toc ，则程序运行完后会显示运行时间。如果你意识到所编的程序运行很慢，需要改进程序，就要注意学习Matlab 变成的规范和技巧。例如在下面的程序中，运用矩阵运算替换了tree.m 中的循环，可以提高运算速度。 改进后的图：

五、程序代码

第一题

pass.m 文件代码

function m=pass(T)

k=50;t=zeros(1,k); %计算k=50辆车，每辆汽车通过路口时间将记录在t 中. tt=solve('-(n-1)*7+(t-N)^2'); %此一元二次方程的2个解中的第1个为正. tt1=11(1); %第n 辆车以加速度通过路口时间tt1(n). ttt=solve('-(n-1)*7+5.5^2+11*(t-n-5.5)'); %第n 辆车以最大限速通过路口时间ttt(n)。 for n=1:k;

t0=5.5+n;x=subs(tt1,n);

if x>t0 %第n 辆车以加速度行驶未达到路口前已达到最大限速. t(n)=subs(ttt,n); %它一定以最大限速通过路口. else t(n)=x; %否则，它以加速度过路口. end end

第二题

图像做旧的程序代码

>> load mask;%读取MATLAB 本身附带的索引图mask ，自动存为ｘ >> mm=X;%mm 数据类型是双精度double ，数值在（0,1）中

>> load face_mos;%读取MATLAB 本身附带的索引图mask ，自动存为X >> nn=X;%nn 数据类型是整数型，数值去0到255

>> mn1=0.5*mm+0.5*double(nn);%都转变成双精度数据，在进行运算 >> mn2=0.8*mm+0.2*double(nn);%取不同的融合比例，得到不同的结果 >> subplot(2,2,1),imshow(mm,[]);

subplot(2,2,2),imshow(nn,[]);subplot(2,2,3),imshow(mn1,[]); subplot(2,2,4),imshow(mn2,[]);%[]作用是拉伸灰度图示图像

第三题

第（1）步：画一个点P （2.5,3）

u=[2.5 3];%用1行2列矩阵表示一个点的坐标

plot(u(1),u(2),'*')%画点，横坐标在前，纵坐标在后，用*号显示

第（2）步：画两个点()3,5.21P

,()6,5.62P >> u=[2.5 3;6.5 6];%用两行两列矩阵表示两个点的坐标

plot(u(:,1),u(:,2),'*')%矩阵第一列u(:,1)表示横坐标，第二列u(:,2)表示纵坐标 axis([2,7,2,7])%给定横坐标范围[2,7]，纵坐标范围[2,7]

第（3）步：画两个点1P

，2P 及其连线上的中点。 u=[2.5 3;6.5 6];

p0=(u(1,:)+u(2,:))/2;%在取中点

u=[u(1,:);p0;u(2,:)];%在原点u中间添入中点坐标

plot(u(:,1),u(:,2),'*-'),axis([2,7,2,7])

第（4）步：

u=[2.5 3;6.5 6];

v=(u(2,:)-u(1,:))/2;%v表示经平移、缩小后的向量OP的终点坐标

A=[cos(pi/6) sin(pi/6);-sin(pi/6) cos(pi/6)]

A =

0.8660 0.5000

-0.5000 0.8660

>> Ip=p0+v*A;%Ip表示经旋转、再平移后的向量P_OP_L的终点坐标

u=[u(1,:);p0;Ip;p0;u(2,:)];

plot(u(:,1),u(:,2),'*-'),axis([2,7,2,7])

第（5）步

branch.m文件代码

function u=branch(v,theta)

p0=(v(2,:)+v(1,:))/2;

v1=(v(2,:)-v(1,:))/2;

A=[cos(theta) sin(theta);-sin(theta) cos(theta)];

lp=p0+v1*A;

u=[v(1,:);p0;p0;lp;p0;v(2,:)];

end

tree.m文件代码

function w=tree(u,theta)

uuu=[];

for j=0;(length(u)/2-1)%ê÷??×üêyê?uDDêy￡¨3¤?è￡?μ?ò?°??￡

u0=u((2*j+1):(2*j+2),:);%òà′?è????óuμ?á?DD113éò???D?μ?￡¨ê÷??￡????óu￡¨￡?

uu=branch(u0,theta);%μ÷ó?éú3éê÷è?μ?×ó3ìDò￡?éú3éò??|ê÷è?

uuu=[uuu;uu];%òàDò?????ùóDê÷è?

end

w=uuu;%μ?μ?±??×??D?μ?ê÷?￡

end

分形树代码

>> clear

u=[0 0;0 1];

subplot(3,3,1),plot(u(:,1),u(:,2)),%逐行逐列画出3行3列9幅图，先将u对应的图放在第1位。

axis([-0.5,0.5,0,1]);%指定图形的坐标范围

for k=1:8

w=tree(u,pi/6);

u=w;

subplot(3,3,k+1),plot(u(:,1),u(:,2)),%将第k阶段生长的树放在第

k+1位

axis([-0.5,0.5,0,1])

end

ans =

ans =

2

ans =

2

ans =

2

ans =

2

ans =

2

ans =

2

ans =

2

改进方法

>> clear

tic

u=[0 0;0 1];

subplot(3,3,1),plot(u(:,1),u(:,2)),axis([-0.5,0.5,0,1]); A=[cos(pi/6) sin(pi/6);-sin(pi/6) cos(pi/6)];

for k=1:8

old=u;

s=length(old)/2-1;

diff=(old(2:2:end,:)-old(1:2:end-1,:))/2;

p1=(old(2:2:end,:)+old(1:2:end-1,:))/2;

p0=p1+diff*A;

u(1:6:6*s+1,:)=old(1:2:end-1,:);

u(2:6:6*s+2,:)=p1;

u(3:6:6*s+3,:)=p1;

u(4:6:6*s+4,:)=p0;

u(5:6:6*s+5,:)=p1;

u(6:6:6*s+6,:)=old(1:2:end,:);

subplot(3,3,k+1),plot(u(:,1)),axis([-0.5,0.5,0,1]) end

toc

Elapsed time is 0.190145 seconds.