数学模型课程设计

数学建模课程设计（程序设计和论文）
题目1.求微分方程的数值解
2.傅立叶级数
3.确定罪犯藏身地点问题
4.曲线拟合与回归分析
5.麦克劳林多项式
6.酒杯的生产
班级/ 学号14140101/2011041401011 学生姓名黄中武
指导教师单锋朱丽梅
沈阳航空航天大学
课程设计任务书
课程名称数学建模实践
院（系）理学院专业信息与计算科学
班级14140101 学号2011041401011 姓名黄中武
课程设计题目1.求微分方程的数值解
2.傅立叶级数
3.确定罪犯藏身地点问题
4.曲线拟合与回归分析
5.麦克劳林多项式
6.酒杯的生产
课程设计时间: 2013年6月17日至2013年7月4日
[要求]
1、学习态度要认真，要积极参与课程设计，锻炼独立思考能力；
2、严格遵守上机时间安排；
3、按照MATLAB编程训练的任务要求来编写程序；
4、根据任务来完成数学建模论文；
5、报告书写格式要求按照沈阳航空航天大学“课程设计报告撰写规范”；
7、报告上交时间：课程设计结时上交报告。

8、严谨抄袭行为。

课程设计的内容及要求：
[内容]
1. 求微分方程的数值解
5
,4,3 1122
11
⎪⎩⎪
⎨⎧=='-='μ='i x q x x q x x x i i
其中初始条件为5,4,3,0)0(,91.46)0(,449.0)0(21====i x x x i 。

其他参数为
5
,4,3 ,/)
12
22*
5
2
22
2
=μ+=μ+=-
⎪⎪⎭
⎫ ⎝⎛+μ=μ∏
=i Y m q Y m q x x k
x x i i i i i m （
,
97.0,69.2,2.2,50,9.360,1026,5.939,2039,67.04322*
5*4*3*2-=-========μm m m k x x x x m 66.11,07.33,69.67,0082.0,26.554325=====Y Y Y Y m
计算在时刻)01.0(,1000,,3,2,1,===h i ih t i 其中 处是函数值。

2. 设函数()f x 是以2π为周期的函数， , 0
()0 , 0x x f x x ππ
-≤<⎧=⎨
≤<⎩，
（1）编写表示函数()f x 的函数M 文件y=fd(x)，并绘图；
（2）求函数()f x 的函数傅立叶级数：)sin cos (2
1
0kx b kx a a
k k k ++
∑
∞
=
（3）对任意的x （x 可以是数组）和n 编写三角多项式)
sin cos (2
)(1
0kx b kx a a
x T k k n
k n ++
=∑
=
的函数M文件y=fly(x,n)
（4）对任意的数组x和n，在同一平面内画出函数(),()
f x T x的图形，并进行比较。

n
3.确定罪犯藏身地点问题
已知罪犯系列犯罪的地点坐标如下：
（1）将角度制的坐标转换为弧度制坐标；
（2）将系列犯罪的地点（13个）的球面坐标转换为平面坐标；并画出这些点的图形（要求带适当宽度的网格，行宽1.25*10^3，列宽h=2.5*10^3）。

（3）求这13个点之间距离的最大值和最小值；分别求这13个点相邻点横坐标和纵坐标之差绝
对值的最大值和最小值；
（4）对所得13个点的每个点，在其周围找16个点（可疑点），其中格宽为d=500。

将这16个点与其余的12个犯罪的地点进行比较，若与某犯罪的地点距离小于等于d ，则去掉该可疑点。

求出所有剩余的点（构成可疑点集合）。

（5）对于d=500:500:60000，重复上述过程，对于不同的d ，求可疑点的个数，画出d 与可疑点个数的图形，并确定d 为何值使，可疑点个数最小。

（6）对上述所得的d ，类似问题4，求可疑点YD 。

（7）对于YD 中每个点（如第j 个点），计算该点到问题2中的13个点(,)i i x y 的距离
,1,2,
,13i s i =，建立该点的衰减函数：13
1
i s j i R e -==∏，求{}j R 的最大值及YD 中对应的点—
—罪犯所在的可能性最大的点，并将该点加在问题2的图形上。

（8）对于d= (21:30)*500，重复上述过程，对于不同的d ，求罪犯所在的可能性最大的点，并将这些点加在问题2的图形上。

4. 曲线拟合与回归分析
(1)以函数M 文件形式编写用1次多项式∑=+=m
j j ij i i x a a y 1
0进行曲线拟合的最小二乘法
通用程序，并根据如下数据计算系数)8,7,6,5,4,3,2,1,0;6,5,4,3,2,1(,==j i a ij .
(2) 再用回归分析法进行线性回归分析
33组
x的值
j
对应的33组
y的值
i
5.（1）编写表示函数()f x 的函数M 文件y=fd(x)，并绘图； （2）求函数)(11ln )(x T n x
x x f n 阶麦克劳林多项式的+-= （3）对任意的x （x 可以是数组）和n 编写多项式()n T x 的函数M 文件y=fly(x,n) 编写对任意固定的n 计算多项式)(x T n 函数值的函数M 文件
（4）对任意的数组x 和n ，在同一平面内画出函数22(),(),[,]33
n f x T x x ∈-的图形，并进行比较。

（5）用()ln 2n T x 计算的近似值，要求误差不超过610-
6. Production of drinking glasses
The main activity of a company in northern France is the production of drinking glasses.
It currently
sells six different types (V1 to V6), that are produced in batches of 1000 glasses, and wishes to plan its production for the next 12 weeks. The batches may be incomplete (fewer than 1000 glasses). The demand in thousands for the 12 coming weeks and for every glass type is given in the following table.
Table : Demands for the planning period (batches of 1000 glasses)
For every glass type the initial stock is known, as well as the required final stock level (in thousands). Per batch of every glass type, the production and storage costs in BC are given, together with the required working time for workers and machines (in hours), and the required storage space (measured in numbers of trays).
The number of working hours of the personnel is limited to 390 hours per week, and the machines have a weekly capacity of 850 hours. Storage space for up to 1000 trays is available. Which quantities of the different glass types need to be produced in every period to minimize the total cost of production and storage?
Table 2: Data for the six glass types
指导教师年月日
负责教师年月日学生签字年月日
沈阳航空航天大学
课程设计成绩评定单
课程名称数学建模实践
院（系）理学院专业信息与计算科学
课程设计题目1.求微分方程的数值解
2.傅立叶级数
3.确定罪犯藏身地点问题
4.曲线拟合与回归分析
5.麦克劳林多项式
6.酒杯的生产
学号2011041401011 姓名黄中武
指导教师评语：
课程设计成绩
指导教师签字
年月日
目录
第一章编程任务 (10)
摘要 (10)
正文 (12)
问题一：求微分方程的数值解 (12)
1.1问题的重述 (12)
1.2问题的分析 (12)
1.3程序框图 (13)
1.4问题的结果 (14)
问题二：傅立叶级数 (14)
2.1问题重述 (14)
2.2问题分析 (14)
2.3程序框图 (15)
2.4问题的结果 (17)
2.4.1问题一的图象 (17)
2.4.2问题二的傅立叶级数证明 (17)
2.4.3问题三的图像 (18)
2.4.4问题四的图像 (19)
问题三：确定罪犯藏身地点 (19)
3.2问题分析 (20)
3.3问题的程序框图 (21)
3.4问题求解 (25)
3.4.1问题一的图象 (25)
3.4.2问题二的结果 (26)
3.4.3问题三的结果 (26)
3.4.4问题四的结果 (27)
3.4.5问题五的图像 (27)
3.4.6问题六的图像(部分) (28)
3.4.7问题七的图像 (29)
3.4.8问题八的图像 (29)
问题四：曲线拟合与回归分析 (30)
4.1问题重述 (30)
4.2问题分析 (32)
4.3程序框图 (32)
4.4问题求解 (33)
4.4.1问题一的求解 (33)
4.4.2问题二的求解 (33)
问题五：麦克劳林多项式展开 (34)
5.1问题重述 (34)
5.2问题分析 (34)
5.3程序框图 (35)
5.4问题的求解 (37)
5.4．1问题一的图像 (37)
5.4．2问题二的麦克劳林多项式 (37)
5.4．3问题三的结果 (38)
5.4．5问题五的结果 (39)
第二章建模任务 (40)
问题六：Production of drinking glasses (40)
6.1问题重述 (40)
6.2问题分析 (40)
6.2问题假设和符号说明 (42)
6.3模型建立 (43)
6.4模型求解 (44)
6.5结果分析与检验 (45)
6.6模型的的优缺点 (46)
参考文献 (48)
第一章编程任务
摘要
对于问题一，求解微分方程组，建立M文件，将题中给出的微分方程组及其所含参数值输入，通过调用函数ode45函数（4-5阶runge-kutta算法，中精度）来实现问题的求解，在命令窗口输入所求范围及初值，从而求出结果。

对于问题二，求函数的傅立叶级数时，首先判断其是否为周期函数，若()f x 不是周期函数，必须先进行周期延拓，但在具体计算过程中，并非用到延拓后的函数()f x 。

在问题（1）中，通过判断x 所在的周期来作出图形；问题（2）在求函数傅立叶变换时，根据()f x 在一个函数周期中的函数表达式求出傅立叶级数的系数；问题（3）将问题（2）所求出的傅立叶级数的系数画出图形；问题（4）将问题（1）和问题（2）所作图形结合在一起，进行比较。

对于问题三，对于问题(1)，我们把度，分，秒独立开来，通过角度之间的换算，最终全部以度作为最后的单位，再将结果除以180/pi ，从而将角度值转换为弧度制。

对于问题(2)，试题中给我们十三组坐标，我们通过建立M 文件，将这十三组坐标列出，通过plot 函数输出。

对于问题(3)，求最小距离和最大距离。

我们首先建立M 文件，将十三组数据输出，通过for 循环，将所有循环出来，再根据两点之间的距离公式算出距离，最后输出最小距离和最大距离。

对于问题(4)，首先建立矩阵存储13给点的坐标，再以中央点为中心构建4*4的小区域，利用循环语句扣除不必要的点，将剩下的点进行距离上的比较，求得记结果。

对于问题(5)，是在上一小题的基础上求解，在此，我们通过距离d 作为循环量，通过不同的d 来画出与可疑点的图形，并在此确定一个d ，是可以点数目最少。

对于问题(6)，对于问题(5)中的d ，按照问题(4)的方法，得到结果。

对于问题(7)，我们按照两点距离公式，计算问题(6)中YD 和问题二中点的距离，并通过建立衰减函数13
1i s j i R e -==∏，求出衰减函数的最大值和
罪犯可能的所在点。

对于问题(8)，在问题(7)的基础上，对不同的d ，以距离d 作为循环量，按照(7)的方法，算出罪犯可能的所在点。

对于问题四，本题考查的是根据多组数据进行最小二乘法曲线拟合，在求取系数时，我们采用《高等代数》上关于最小二乘法的求解方法，通过矩阵之间的运算求出系数。

对于问题五，根据()f x 麦克劳林的推导公式，我们可以用matlab 求出它的n 次麦克劳林展式并将其以表达式()n T x 的形式输出出来。

然后，编写M 文件求对于给定的次数n 和x 求出()n T x 的值。

最后对于第三小问，用plot 命令将函数()f x 的值与其的麦克劳林展式()n T x 在区间22
[,]33
上进行比较。

关键词：微分方程数值解；傅立叶级数；角度与弧度；最小二乘法；
正文
问题一：求微分方程的数值解
1.1问题的重述
5,4,3 112211⎪⎩⎪⎨⎧=='-='μ='i x q x x q x x x i i
其中初始条件为5,4,3,0)0(,91.46)0(,449.0)0(21====i x x x i 。

其他参数为
5
,4,3 ,/)1222*5222
2=μ+=μ+=-⎪⎪⎭⎫ ⎝⎛+μ=μ∏=i Y m q Y m q x x k x x i i i i i m （
,
97.0,69.2,2.2,50,9.360,1026,5.939,2039,67.04322*5*4*3*2-=-========μm m m k x x x x m 66.11,07.33,69.67,0082.0,26.554325=====Y Y Y Y m
计算在时刻)01.0(,1000,,3,2,1,===h i ih t i 其中 处是函数值。

1.2问题的分析
求解微分方程组，建立M 文件，将题中给出的微分方程组及其所含参数值输入，通过调用函数ode45函数（4-5阶runge-kutta 算法，中精度）来实现问题的求解，在命令窗口输入所求范围及初值，从而求出结果。

H
1.3程序框图
图1.3 问题一的程序框图
1.4问题的结果
图1.4 问题一的结果
问题二：傅立叶级数
2.1问题重述
设函数()f x 是以2π为周期的函数， , 0
()0 , 0x x f x x ππ-≤<⎧=⎨≤<⎩
，
（1）编写表示函数()f x 的函数M 文件y=fd(x)，其中x 可以是任意实数，也可以是任意数组，并在任意区间[,]a b 上绘出函数()f x 图形；
（2）推倒函数()f x 的函数傅立叶级数：)sin cos (2
1
0kx b kx a a
k k k ++
∑
∞
=
(3)对任意的x （x 可以是数组）和n 编写三角多项式
)sin cos (2
)(1
0kx b kx a a
x T k k n
k n ++
=∑
=的函数M 文件y=fly(x,n)
（4）对任意的数组x 和n ，在任意区间[,]a b 上在同一平面内画出函数(),()n f x T x 的图形，并进行比较。

2.2问题分析
对于问题一，通过建立M 文件，输入自变量的范围和函数的表达式，通过plot 函数画图。

对于问题二，对于傅立叶基数的证明，可参照《数学分析》上的证明。

对于问题三，该函数为周期函数，其图像为周期图像，首先建立M 文件，输入自变量的范围，通过for 语句的循环，来实现傅立叶级数程序的运行，最后通过plot 函数，来实现图形的绘出。

对于问题四，建立M 文件，在M 文件中分别写上两个函数图像的运行程序，通过hold on 来实现统一窗口绘出两个函数的图像。

2.3程序框图
图2.3-1 问题一y=fd(x)程序框图
图2.3-2 问题三y=fly(x,n)程序框图
2.4问题的结果
2.4.1问题一的图象
图2.3.1 fd(x)的函数图像
开始
调用函数fd （x ，n ）
输入x ，n
调用函数
调用plot 函数
结束
图2.3-3 同时绘制()f x 与其傅里叶级数图像的程序框图
2.4.2问题二的傅立叶级数证明
求函数()f x 的函数傅立叶级数：)sin cos (2
1
0kx b kx a a
k k k ++
∑
∞
=
证明：函数展开成傅里叶级数，设)(x f 是定义在)π,π(-上的有界周期函数，则能展开成三角级数
∑∞
=++=1
)
sin cos (2)(k k k kx b kx a a x f .)1(0a 求x
kx b kx a x a x x f k k k d ])sin cos ([d 2d )(ππ
1
π
π0
π
π⎰∑⎰⎰
-∞=--++=x kx b x kx a x a k k k k d sin d cos d 2ππ
1
ππ1π
π0
⎰∑⎰∑⎰-∞=-∞=-++=,22
π⋅=
a .
d )(π1π
π0⎰-=x x f a 则]
d cos sin d cos cos [π1
⎰⎰∑--∞
=++π
π
π
x nx kx b x nx kx a k k k ⎰-=π
π
x
nx a n d cos 2,
πn a =⎰-
=
π
ππ
x nx x f a n d cos )(1
则).
,3,2,1( =n ⎰⎰-
-
=π
ππ
π
x
nx a x nx x f d sin 2
d sin )(0
]
d sin sin d sin cos [1
⎰⎰∑--∞
=++π
π
ππ
x nx kx b x nx kx a k k k .
)3(n b 求.
)
2(n a 求⎰⎰-
-
=
π
ππ
π
x
nx a x nx x f d cos 2
d cos )(0
综上求得结果为：
1
021(1),(1cos ),2n n n a a n b n n
π
ππ+-=-=-=
2.4.3问题三的图像
图2.3.3 y=fly(x,n)函数图像
2.4.4问题四的图像
,
πn b =).
,3,2,1( =n
图2.3.4 20阶傅里叶图像与原图像比较
问题三：确定罪犯藏身地点
3.1问题重述
已知罪犯系列犯罪的地点坐标如下：
（1）将角度制的坐标转换为弧度制坐标；
（2）将系列犯罪的地点（13个）的球面坐标转换为平面坐标；并画出这些点的图形（要求带适当宽度的网格，行宽1.25*10^3，列宽h=2.5*10^3）。

（3）求这13个点之间距离的最大值和最小值；分别求这13个点相邻点横坐标和纵坐标之差绝对值的最大值和最小值；
（4）对所得13个点的每个点，在其周围找16个点（可疑点），其中格宽为d=500。

将这16个点与其余的12个犯罪的地点进行比较，若与某犯罪的地点距离小于等于d ，则去掉该可疑点。

求出所有剩余的点（构成可疑点集合）。

（5）对于d=500:500:60000，重复上述过程，对于不同的d ，求可疑点的个数，画出d 与可疑点个数的图形，并确定d 为何值使，可疑点个数最小。

（6）对上述所得的d ，类似问题4，求可疑点YD 。

（7）对于YD 中每个点（如第j 个点），计算该点到问题2中的13个点(,)i i x y 的距离
,1,2,
,13i s i =，建立该点的衰减函数：13
1
i s j i R e -==∏，求{}j R 的最大值及YD 中对应的点—
—罪犯所在的可能性最大的点，并将该点加在问题2的图形上。

（8）对于d= (21:30)*500，重复上述过程，对于不同的d ，求罪犯所在的可能性最大的点，并将这些点加在问题2的图形上。

3.2问题分析
对于问题(1)，我们把度，分，秒独立开来，通过角度之间的换算，最终全部以度作为最后的单位，再将结果除以180/pi ，从而将角度值转换为弧度制。

对于问题(2)，试题中给我们十三组坐标，我们通过建立M 文件，将这十三组坐标列出，通过plot 函数输出。

对于问题(3)，求最小距离和最大距离。

我们首先建立M 文件，将十三组数据输出，通过for 循环，将所有循环出来，再根据两点之间的距离公式算出距离，最后输出最小距离和最大距离。

对于问题(4)，首先建立矩阵存储13给点的坐标，再以中央点为中心构建4*4的小区
域，利用循环语句扣除不必要的点，将剩下的点进行距离上的比较，求得记结果。

对于问题(5)，是在上一小题的基础上求解，在此，我们通过距离d 作为循环量，通过不同的d 来画出与可疑点的图形，并在此确定一个d ，是可以点数目最少。

对于问题(6)，对于问题(5)中的d ，按照问题(4)的方法，得到结果。

对于问题(7)，我们按照两点距离公式，计算问题(6)中YD 和问题二中点的距离，并通过建立衰减函数13
1
i s j i R e -==
∏
，求出衰减函数的最大值和罪犯可能的所在点。

对于问题(8)，在问题(7)的基础上，对不同的d ，以距离d 作为循环量，按照(7)的方法，算出罪犯可能的所在点。

3.3问题的程序框图
图3.3-1 问题一坐标转换流程图
图3.3-2 问题二的绘制网格的流程图
图3.3-3 问题三的最值流
图3.3-4 问题四的距离上的比较流程图
图3.3-5 求解总个数最小时所对应d的流程图
图3.3-6 求可疑点流程图
图3.3-7 画图求点流程
附加说明：由于问题(5)和问题(6)算法一样，所以流程图可以共用一个
3.4问题求解
3.4.1问题一的图象
图3.4.1 问题一的结果
3.4.2问题二的结果
问题二的图象：
图3.4.2 问题二的图像3.4.3问题三的结果
图3.4.3 问题三的图像
3.4.4问题四的结果
图3.4.4 问题四的结果3.4.5问题五的图像
图3.4.5 问题五的图像
3.4.6问题六的图像(部分)
图3.4.6 问题六的结果
3.4.7问题七的图像
图3.4.7 平面坐标上的可疑点
3.4.8问题八的图像
图3.4.8 罪犯所在的可能性最大的点
问题四：曲线拟合与回归分析
4.1问题重述
(1)以函数M 文件形式编写用1次多项式∑=+=m
j j ij i i x a a y 10进行曲线拟合的最小二乘法
通用程序，并根据如下数据计算系数)8,7,6,5,4,3,2,1,0;6,5,4,3,2,1(,==j i a ij .
(2)再用回归分析法进行线性回归分析。

33组j x 的值
对应的33组
y的值
i
4.2问题分析
本题考查的是根据多组数据进行最小二乘法曲线拟合，在求取系数时，我们采用《高等代数》上关于最小二乘法的求解方法，通过矩阵之间的运算求出系数。

4.3程序框图
图4.3 最小二乘法曲线拟合的流程图
4.4问题求解
4.4.1问题一的求解
图4.4.1 问题一的结果4.4.2问题二的求解
图4.4.2 问题二的结果
问题五：麦克劳林多项式展开
5.1问题重述
（1）编写表示函数()f x 的函数M 文件y=fd(x)，并绘图； （2）求函数)(11ln
)(x T n x
x
x f n 阶麦克劳林多项式的+-= （3）对任意的x （x 可以是数组）和n 编写多项式()n T x 的函数M 文件y=fly(x,n) 编写对任意固定的n 计算多项式)(x T n 函数值的函数M 文件
（4）对任意的数组x 和n ，在同一平面内画出函数22
(),(),[,]33
n f x T x x ∈-的图形，并
进行比较。

（5）用()ln 2n T x 计算的近似值，要求误差不超过610-。

5.2问题分析
针对本题，我们需要将x
x
x f +-=11ln
)(转化为)1ln()1ln()(x x x f +--=，在此，我们对)1ln()(x x A -=和)1ln()(x x B +=分别进行麦克劳林展开为 )(1x T n 和)(2x T n ，则
)()()(x B x A x f -=展开后可表示为)()()(21x T x T x T n n n -=，编写M 文件求出其表达式
形式即可。

对于得到的表达式结果用Matlab 利用字符表达式带入值很容易求出结果。

在第三问中，输入
x 的值将其带入)(x f 和)(x T n 中，并求出它们的差值)(x E n 。

为了方便画图，将程序做一修改，使其可以输入一个x 向量，可以得到一个结果向量。

然后对其在区间
]32
,32[ 画图进行比较。

5.3程序框图
图5.3.1 问题一的函数流程图。