matlab程序大作业

中国石油大学（北京）2017 — 2018 学年第一学期《Matlab语言》期末程序设计报告

班级：

姓名：

学号：

一.(10分) 有砂岩含气和含水情况下的速度频散曲线如图1a、1b所示(数据请到网络课程下

载)，请用Logistic函数对两条曲线进行拟合。要求: 1)借用超定线性方程求解实现拟合，避免使用fit、polyfit等函数；2）给出拟合所采用Logistic函数形式；3) 给出针对两种情形进行拟合得到的参数值；4）绘出拟合结果；5）代码附于附录1。

（a）（b）

图1.某砂岩含气（a）和含水（b）情况下的速度频散曲线。

答: 2）

3）由下图得来k = 0.07

5

4）油 水

frequency

v g a s

frequency

v w a t e r

5)见附录一

二.（30分）多位学者研究表明，地震反演的分辨率极限与震源子波的波峰波谷间隔以及拐点间隔大小紧密相关，而波峰波谷以及拐点位置又是峰值频率f 的函数。现已知零相位

雷克子波的表达式为2

)(2)

)(21()(ft e ft t w ππ--=， a) 求该子波的极值（即波峰与波谷）位置t 1，给出t 1的解析表达式，并计算当频率f = 30 Hz 时t 1的取值；

b) 求该子波的拐点位置t 2，给出t 2的解析表达式，并计算当频率f = 30 Hz 时t 2的取值；

请将代码附于附录2。 答：a ）t1 =

/（2*pi*f ），0，-/（2*pi*f ）

当f = 30时，t1 = 0.0130，0，-0.0130

b ）t2 =

/（2*pi*f ），

-

/（2*pi*f ），

/（2*pi*f ），

-

/（2*pi*f ）；

当f = 30时，t2 = 0.0175，-0.0175，0.0056，-0.0056 c ）程序见附录2

三. (60分)已知乐透型彩票单注金额为2元，一般的销售规则为：低项奖数额固定，高项

奖按比例分配，各高项奖额的计算方法为：

[当期销售总额*总奖金比例 – 低项奖总额] * 单项奖比例。

其中总奖金比例一般为0.5，一、二、三等奖为高项奖，四、五、六、七为低项奖。现准备在某地区投入m/n （即从n 个号码中任意挑选m 个号码）型乐透彩票，已知该地区平均收入和消费水平的相关因子α=630589，并且彩民对奖金x 的满意度为：

2

)(e

1)(α

x x u --=

要寻求该地区的最佳销售策略。该问题抽象为在满足一定约束条件下使目标函数达到最

大的非线性规划问题。 目标函数为：

∑==

7

1

)(i i

i

x u p F

约束条件为：

3,2,1,)1(7

4

=-=

∑=j p x p x j

j

i i

i j r

（3.1）

1321=++r r r

8.05.01≤≤r

61510*510*6≤≤x

6,...,2,1,1

=≤≤

+i b x x a i i i

i 6,...,2,1,1=<+i p p i i

为整数,75m m ≤≤

为整数，6029n n ≤≤

0,0≥>j j x r

其中，r j 为j 等奖占高项奖的比例，j =1,2,3;x i 为i 等奖奖金额；p i 为中i 等奖的概率，i =1，2，…，7；a 1=10,b 1=233, a 2=4, b 2=544, a 3=3, b 3=17, a 4=4, b 4=20, a 5=2, b 5=20, a 6=2, b 6=20。

已知计算m/n 型乐透彩票中1至7等奖概率的代码如下：

function p = cal_pro(m, n) % 计算各奖项的概率

p = zeros(7, 1);

c = comb(n, m);

p(1)=1/c;

p(2)=comb(m,m-1)./c;

p(3)=comb(m,m-1) * comb(n-m-1,1) / c;

p(4)=comb(m,m-2) * comb(n-m-1,1) / c;

p(5)=comb(m,m-2) * comb(n-m-1,2) / c;

p(6)=comb(m,m-3) * comb(n-m-1,2) / c;

p(7)=comb(m,m-3) * comb(n-m-1,3) / c;

end

function c = comb(n, m)

% 从n个数中取出m个数的组合

c = factorial(n) / factorial(n - m) / factorial(m);

end

请用Matlab编程求解该地区销售乐透型彩票的最佳方案，即求m、n、r1、r2、r3、x4、x5、x6、x7的最优解。要求：1）请给出m、n、r1、r2、r3、x4、x5、x6、x7的取值；2）代码（不包含给出的cal_pro和comb的代码）附于附录3。

答：(1).m = 43,n = 6,r1 = 0.5669,r2 = 0.2981,r3 = 0.1350,x4 = 317.9144,x5 = 40.2352,x6 = 11.6676,x7 = 5.7963

(2)见附录三

附录1

clc;

clear all;

close all;

load v_f_gas;%载入油文件

v_gas = v;

figure;

plot(f,v_gas,'b-');hold on;

xlabel('frequency');

ylabel('v_gas');

k = [-0.5:0.01:0.5];