数学建模课后习题答案

方程及方程组的求解
路灯照明问题。

在一条20m 宽的道路两侧, 分别安装了一只2kw 和一只3kw 的路灯, 它们离地面的高度分别为5m 和6m 。

在漆黑的夜晚, 当两只路灯开启时 （1）两只路灯连线的路面上最暗的点和最亮的点在哪里？ （2）如果3kw 的路灯的高度可以在3m 到9m 之间变化, 如何路面上最暗点的亮度最大？ （3）如果两只路灯的高度均可以在3m 到9m 之间变化, 结果又如何？
解:
根据题意, 建立如图模型
P1=2kw P2=3kw S=20m 照度计算公式:
2
sin r p k I α
= （k 为照度系数, 可取为1；
P 为路灯的功率）
（1）设Q(x,0)点为两盏路灯连线上的任意一点, 则两盏路灯在Q 点的照度分别为
21111sin R p k I α= 222
22
sin R p k I α=
221
21
x h R += 1
11sin R h =
α
2
22
2
2)(x s h R -+= 2
2
2sin R h =
α
Q 点的照度:
3
23
23
222
2
23
221
11))20(36(18)
25(10)
)((()
(()(x x x s h h P x h h P x I -++
+=
-++
+=
要求最暗点和最亮点, 即为求函数I(x)的最大值和最小值, 所以应先求出函数的极值点
5
25
25
222
225
221
11'
))20(36()20(54)
25(30)
)(()(3)
(3)(x x x x x s h x s h P x h x h P x I -+-+
+-=
-+-+
+-=
利用MATLAB 求得0)('=x I 时x 的值
代码:
s=solve('(-30*x)/((25+x^2)^(5/2))+(54*(20-x))/((36+(20-x)^2)^(5/2))'); s1=vpa(s,8); s1
运行结果: s1 =
19.97669581 9.338299136 8.538304309-11.61579012*i .2848997038e-1
综上, x=9.33m 时, 为最暗点；x=19.97m 时, 为最亮点。

路灯2的高度可以变化时, Q 点的照度为关于x 和h2的二元函数:
3
222
2
3
23
222
223
221
112)
)20((3)
25(10)
)(()
(),(x h h x x s h h P x h h P h x I -++
+=
-++
+=
与（1）同理, 求出函数I(x,h2）的极值即为最暗点和最亮点
0)
)((3))((52222
2
2322222=-+--+=∂∂x s h h P x s h P h I 利用matlab 求得x:
solve('3/((h^2+(20-x)^2)^(3/2))-3*(3*h^2)/((h^2+(20-x)^2)^(5/2))=0') ans =
20+2^(1/2)*h 20-2^(1/2)*h
即x1=20+2^(1/2)*h （舍去） x2=20-2^(1/2)*h
0)
)20(()20(9)25()220(30-))(()(3)(35222252522222522111=-+-++-=-+-++-=∂∂x h x h x h x s h x s h P x h x h P x I
利用matlab 求解h 2
solve('-30*(20-2^(1/2)*h)/((25+(20-2^(1/2)*h)^2)^(5/2))+9*h*(20-(20-2^(1/2)*h))/((h^2+(20-(20-2^(1/2)*h))^2)^(5/2))=0') ans =
7.4223928896768612557104509932965 14.120774098526835657369742179215 因为h 在3~9之间, 所以h2=7.42239m 再利用matlab 求解x 和亮度I
算法: h=7.42239; x=20-2^(1/2)*h
I=10/((25+x^2)^(3/2))+(3*h)/((h^2+(20-x)^2)^(3/2)) 结果: x =
9.5032 I =
0.0186
综上, x=9.5032 , h2=7.42239时, 最暗点的亮度最大, 为0.0186w 。

两盏路灯的高度均可以变化时, I 为关于x,h1,h2的三元函数, 用同样的方法求解
3
222
2
23
221
1121)
)(()
(),,(x s h h P x h h P h h x I -++
+=
0)
(3)(52212
1
1322111=+-+=∂∂x h h P x h P h I
)
)20((9))20((3))((3))((52222
2
322252222
22322222=-+--+=-+--+=∂∂x h h x h x s h h P x s h P h I
0)
)20(()20(9)(6))(()(3)(35222252211522222522111=-+-++-=-+-++-=∂∂x h x h x h x h x s h x s h P x h x h P x I
x
h
2
1
1
=
)
20
(
2
1
2
x
h-
=
2
5
2
2
1
1
2
5
2
2
2
2
)
(
2
]
)
20
(
[
)
20
(
3
x
h
x
h
x
h
x
h
+
=
-
+
-
2
5
2
2
2
2
5
2
2
2
)
2
1
(
2
2
]
)
20
(
)
20
(
2
1
[(
)
20
(
2
3
x
x
x
x
x
x
+
=
-
+
-
-
=
3
33
2
)
20
(
1
x
x
=
-
利用matlab求解x, h1, h2的值:
算法: solve('1/((20-x)^3)=2/(3*(x^3))');
s1=vpa(s,6);
a=(1/sqrt(2))*s1;
a1=double(a);
b=(1/sqrt(2))*(20-s1);
b1=double(b);
a1,b1,s1
结果:
a1 =
6.5940
5.1883 +12.0274i
5.1883 -12.0274i
b1 =
7.5482
8.9538 -12.0274i
8.9538 +12.0274i
s1 =
9.32530
7.33738+17.0093*i
7.33738-17.0093*i
综上, h1 =6.5940, h2=7.5482 , x=9.32530时, 最暗点的亮度最大
数据插值
山区地貌: 在某山区测得一些地点的高程如下表3.8。

平面区域为
（1200<=x<=4000,1200<=y<=3600) 试作出该山区的地貌图和等高线图, 并对几种插值方法进行比较。

表3.8 某山区高程表
利用matlab编程代码如下:
x=1200:400:4000;
y=1200:400:3600;
[xi,yi]=meshgrid(1200:4000,1200:3600);
z=[1130 1250 1280 1230 1040 900 500 700;
1320 1450 1420 1400 1300 700 900 850;
1390 1500 1500 1400 900 1100 1060 950;
1500 1200 1100 1350 1450 1200 1150 1010;
1500 1200 1100 1550 1600 1550 1380 1070;
1500 1550 1600 1550 1600 1600 1600 1550;
1480 1500 1550 1510 1430 1300 1200 980];
线性插值法
zi=interp2(x,y,z,xi,yi,'linear');
mesh(xi,yi,zi)
title('线性插值法')
xlabel('x');
ylabel('y');
zlabel('z');
C=contourf(xi,yi,zi);
clabel(C);
title('等高线图')
最邻近插值法
zi=interp2(x,y,z,xi,yi,'nearest');
mesh(xi,yi,zi)
title('最邻近插值法')
xlabel('x');
ylabel('y');
zlabel('z');
C=contourf(xi,yi,zi);
clabel(C);
title('等高线图')
立方插值法
zi=interp2(x,y,z,xi,yi,'cubic' );
mesh(xi,yi,zi)
title(' 立方插值法')
xlabel('x');
ylabel('y');
zlabel('z');
C=contourf(xi,yi,zi);
clabel(C);
title('等高线图')
三次样条插值法
zi=interp2(x,y,z,xi,yi,'spline' );
mesh(xi,yi,zi)
title(' 三次样条插值法')
xlabel('x');
ylabel('y');
zlabel('z');
C=contourf(xi,yi,zi);
clabel(C);
title('等高线图')
四种差值方法在
运算时间光滑程度运算时间和光滑
程度上有一定的
差异, 如下表所
示
类别
差值方法
最邻近插值法快差
线性插值法稍长稍好三次样条插值法最长最好
曲线拟合
方法一
利用1stOpt快速拟合公式搜索可得到公式为:
y = p1+p2*x+p3/x+p4*x^2+p5/x^2+p6*x^3+p7/x^3+p8*x^4+p9/x^4+p10*x^5 p1=18382690.6773727
p2=-4152096.11663013
p3=-51037385.3263795
p4=592195.144413008
p5=84947107.1889704
p6=-51716.5130172659
p7=-75932896.2582835
p8=2521.12152863706
p9=27252247.5649699
p10=-52.482670759974
Matlab代码如下
p1=10802.6249167589;
p2=-20010.6348923663;
p3=19400.634311511;
p4=-10100.4704562703;
p5=2958.58084727337;
p6=-461.436321152701;
p7=21.9610124897453;
p8=4.50124440221874;
p9=-0.851576261728162;
p10=0.0575464303972622;
p11=-0.00144545223415816;
x=4.5;
y=p1+p2*x+p3*x^2+p4*x^3+p5*x^4+p6*x^5+p7*x^6+p8*x^7+p9*x^8+p10*x^9 +p11*x^10
运行结果:
y =
921.1616
方法二
利用matlab拟合
x=[1 2 3 4 5 6 7 8 9 10];
y=[2615 1943 1494 1087 765 538 484 290 226 204];
plot(x,y,'*')
观察图形可知, 曲线近似于指数函数
设 , 取对数得
记a a k a y m ln ,,ln 21=-==,则21a x a m += 利用matlab 算出a1,a2的值
x=[1 2 3 4 5 6 7 8 9 10];
y=[2615 1943 1494 1087 765 538 484 290 226 204]; m=log(y);
aa=polyfit(x,m,1)
结果
aa =
-0.2969 8.1591
则8.1591-0.2969x
m += , 即1591.8-0.2969x ln +=y 所以1591.8-0.2969x +=e y
利用matlab 算出x=4.5时y 的值
x=4.5;
y=exp(-0.2969*x+8.1591)
结果
y =
918.7830。