中国石油大学计算实习报告

2017-2018-2学期《计算实习》实习报告
1、实习总结
我大二曾修读过崔老师的《 Matlab 程序设计》，但是没想到这次计算实习用起来已经生疏了不少，时常需要参考课本。

此外，这门课的难度也颇高，许多问题都没有固定的解答套路，所以经常需要我和同学交流讨论想法，使自己的算法进一步完善。

好在同桌和好朋友也不吝惜自己的想法，时常使我茅塞顿开。

经过本次实习，为了更好地总结学习到的各种编程方法与思想，写了这份实习报告。

总得来说我的收获分为 3 方面：
1.应用是理论的反馈
实习题目2.2中证明两个向量组等价，我的第一反应就是"存在一组不全为 0的” 并且将系数解出来的方法，然而这种固定得思维模式却不利于解决后面的问题，崔老师通过简单的移项就证了出来，让我有点打开了新世界的感觉，豁然开朗，也提醒我不要固化思维。

实习题目 2.3 中，一开始的求导就有点让我寸步难行，这是由于我对链式求导法掌握程度不够熟练的原因，也就是数学基础不够扎实。

实习题目 2.8 中，一开始我也很大意，也想直接利用牛顿 -莱布尼兹公式求积分，但却忽视了 y 和 x 潜在的某种函数关系，所以我们不能直接积分。

我又发现求积节点等距分布，所以我很自然得考虑利用数值分析中的复化Simpson公式，接着又想到可以利用梯形公式进一步优化自己的算法。

再比如，我们要加强类似题目之间的相关性、串联性，比如题 2.9 和 2.11 都是需要用到四阶龙格库塔方法，如果把题 2.11 原理弄明白了，那就能很快得到 2.9 的思路和算法。

诸如此类的感想我觉得都可以归结到应用能力对基础理论知识的反映的这一类问题上。

我在实习的时候时常因为基础数学理论知识不够扎实，导致反馈到应用计算上面时就经常会出现方法上的问题。

这就需要我更加注重基础数学知识的打磨和积累，在编程的时候需要更加注重对理论知识的理解。

2.加强对理论本质的理解
Matlab 自带的工具包虽然能让我们更加快速高效得编写程序，但也容易让我们养成依赖，而且如果我们过于依赖现成的工具，而不了解理论本质，就容易陷入固化的思维模式中。

正如最近的中美贸易战和“中兴事件”都是由于我们没有掌握“核心”而被人掐住了要害，受制于人。

所以在本次计算实习的时候，崔老师许多题目不再让我们直接使用一些简单的函数命令，比如实习题目 2.1 等我们本可以直接利用“ ”这个运算符，但我们从本质出发，编写了列主元高斯消去法，这有利于进一步加深了对求解方程组的相关知识的理解。

3.团队合作与虚心请教
由于我的数学基础知识不如刘劭胡汇丰扎实，所以需要我常常课后和同学们探讨思路和算法，目的是集思广益，寻找各方的优缺点，进而设计一种可编写的、在我本人能力范围内的算法。

此外，这次计算实习虽然难度比较大，但对我来说是在大学里最后一次能很好提升自己的机会（因为社会上更加注重团队合作，很多人也不乐意帮你解答困惑）。

以前我连很多函数命令都用得不熟悉，在 7月 5日课程结束后我花费了大量时间，潜心坐在电脑前思考，并虚心向许多同学请教，不断调试，终于完成了这次实习。

当然这次实习不只是为了这份实习报告，它也磨练我的心智，更重要得是把这些实习内容内化成了自己的编程能力。

未来我需要积极培养创新型的思维模式，努力成为复合型人才。

2、实习过程
2.1实习题目1 题目
用Crank - Nicholson 格式求解一维热传导方程的初边值问题
2
汀一 ::T
一 ；x 2
求解过程或算法：
首先参考《数值分析》 P67的方法，编写了列主元高斯消去法的
matlab 程序：
function [X] = Gauss(B)
[u,v]=size(B); %求出输入矩阵 B 的维度 m=u;
X=zeros(m,1) %类似一个初始化，建立一个 m 行1列的0矩阵
for n=1:1:m-1 %一个循环嵌套，步长为
1，从1到m-1行(因为最后一行 m 行对角线上元素不取)
u=B(n:m,n); %建立一个矩阵n 行到m 行，位置第n 列的一个(m-n )行1列的矩阵 [c,i]=max(u); %求里面的最大元素的
位置
if c==0 %如果列主元为 0，贝U detA=0
error( 'error'); else
if i+n-1~=n
%如果不在对角线上
t=B(n,:); %取一个中间变量，最大元所在 n 行的一整行元素
B(n,:)=B(i+n-1,:);
B(i+n-1,:)=t; %把列主元所在的第 n 行和第i+n-1行整行对换 end end
for j=n+1:1:m
B(j,:)=B(j,:)-B(j,n)./B(n,n)*B(n,:); %对角元以下消元计算，变成一个上三角矩阵
0 ：： x :: 1,t 0
求出t = 0.2时刻的温度分布。

<T (x,0)=sinux (O c x c l ) T 0,t = T 1,t = 0
T n
n n
__ j+一
2T
j +T
j 」 =
2
(
h 2 T ；i 1
_
2T j
n1
j 1
、
h )
也可写为
n 1
n
1
n n n n ：1 n+1 n+1
Tj =Tj +J((Tj *—2Tj +Tj _J+(Tj 卅-2Tj
))
其中
=亍001=1
h 0.01
(1-1 )
(1-2)
(1-3)
end end
A=B(:,1:m); %1到m 列的数据
d=B(:,m+1); %增广矩阵最后一列，即 b for k=m:-1:1 % 倒序
z=0;
for l=1:1:m-k
z=z+A(k,k+l)*X(k+l); %求接上三角方程组
end
X(k)=(d(k)-z)/A(k,k); % 回代求解
end end
再对Crank-Nicholson 格式对n 和n+1移项化简，得到
(1 +呻1
・話皆-窪禺=U -呵+轨+1喑
所以我们可以得到等式两侧的矩阵形式:
根据上面的三对角阵和第 n 层的T 值，运用高斯消元法求得第 n+1层的T 值，重复此 过程可求的一段时间内的 T 值，进而可画出t 时刻的数值解的图像。

再根据解析式求得t=0.2 时刻的精确解。

于是我们可以编写
t=0.2时的CN 格式程序：
function [u]=CN(A,T) %CN 格式，A 为网格比，T 为对时间的划分 x=0:0.1:1; m=l ength(x);
z=zeros (T+1,m); % 温度矩阵 u=si n( pi*x);z(1,:)=u;
C=(1+A)*eye(m-2)+diag((-A/2)*ones(1,m-3),1)+diag((-A/2)*ones(1,m-3),-1); %我给出的三对角阵写法 对角线上元素+对角线上1行+对角线-1行(diag 函数) for n=1:1:T
u(1)=0;u(m)=0;
b=(1-A)*u(2:m-1)+A/2*u(3:m)+A/2*u(1:m-2); %将上面等号右侧的矩阵写成这样的格式
Z=[C,b'];
u(2:m-1)=Gauss(Z); %用自己编写的列主元高斯消元法求解 z(n+1,:)=u; %对温度赋值 end
/1 + A -钛
"・
■an
0 Q
Q \
% £■»
1-A
ELV
■
E
■U
%
-X
3
1-V
x=0:0.1:1;X=0:0.05:1;t=0:0.1*0.1*A:0.2;
U=exp(-piA2*0.2)*sin(pi*X); %t=0.2 时精确解
plot(x,u,'black' ,X,U：red')%画精确解与数值解的二维比较图
legend('入=1':精确解',4)
title('Crank-Nicolson 格式二维图像');
figure
surf(x,t,z) %画数值解的三维图
title('Crank-Nicolson 格式三维图像');
xlabel('间距x');
ylabel('时间t');
zlabel('温度z');
程序结果：
在comma nd wi ndow里运行CN(1,20)，我们可以得到如下两张结果图:
图1 图2
结论：由图1可以看出，在网格比/二：［时，得出的数值解与精确解非常近似，这是因为C-N 格式是无条件稳定的。

2.2实习题目2
题目设v是非零n维向量，A是n xn可逆矩阵，与 A相关的向量组为
v, A V,A2V,III, A k4v (2-1)
线性无关，记 K k=span{ v, Av, A2v, ..., A k-1v}。

以k = 4为例，如果按如下的计算过程：
①将v单位化：令V1 = v/||v||,然后用v1替换原向量组中的 v,得
2 3
v1, Avi, A v-!, A v1(2-2)
该向量组显然与原向量组等价。

②利用Gram - Schmidt正交化方法，从 Av1中减去它在v1上的投影向量，得
v2 = Aw - Av1 ,vi v(2-3)
单位化后得V 2=V 2/||V 2|。

用V 替换向量组（2-2）中的AV 1,得:
V i , V 2, Av,, A 2
V 2
（1） 请你证明：向量组（I ）和（II ）等价，即它们可以相互线性表示。

（2） 请你按照①、②的逻辑过程，继续对向量组进行执行计算、替换操作，写出计算 过程，最终得到如下的规范化正交向量组（标准正交基）
V i , V 2,V 3, V 4
h
i,j 二
AV j ,
V i
, h j i,j = V j i
（4） 请你计算（归纳）出以上算法的矩阵形式（提示：请参考利用 Gram - Schmidt 正
交化方法对矩阵进行 QR 分解的计算过程）。

（5） 请讨论以上的计算过程与 Gram - Schmidt 正交化方法之间的关系。

求解过程或算法：
1）我们先证（n ）可以由（I ）线性表示出来： 我们对n.「移项化简并带入「=‘；川，得到下式
血國卜（血曲國E
记为冷衣%，錶尿
两边同乘一，得到: - ：_ - : 同理可得 ，得到二一 _
_
•••（n ）可以由（I ）线性表
示 再证（I ）可以由（n ）线性表示：
类似的，我们对 试厂理％…软蚁褥酗•.移项化简得到并带入 一表达式即
航讦叫1冊11 +
（加”叫）耳
记为• 罰&品能
对上式两边同乘二，得到_ - -一 同理可得，两边同乘护，得到.
......
• （I ）可以由（n 熾性表示
综上得：向量组（I ）和向量组（n ）可以互相线性表出，即是等价的。

证毕
2） 利用Gram-Schmidt 正交化方法，从空％:中减去它在馬「：底上的投影向量，得
（3）请你将以上计算过程归纳为以
k 为循环变量的一个算法。

提示：引进新的记号:
(2-4)
(2-5)
(2-6)
. —
单位化后得&:;-:川〔；|。

用•一替换向量组（I ）中 的血雄，得：一.殴飞麹.
利用Gram-Schmidt 正交 化方法，从」i 中减去它在说| ：：险 險上的投影向量，得
二一
. 二］:.:,单位化后得
用_替换向
量组（I ）中的姒；,得： 3）算法：
1
.开始：._=』|忌||
2•迭代
for j=1 . k-1 do
hjj =巧=1,2 ....... j
i
町+1=阿_2如円
i
End 4）将第二问的结果写成矩阵形式，结合第三问的结果：
■
观察算法中式子的'
以及占曲隔;.得到
矽为上下两块©）心釉吩别梯怖和最后-彳亍
叫7 陋,H k 为k-1阶方阵，且胪的前k-1个元素都为零
^k-i>fc-v
•••归纳出的以上算法的矩阵形式为 厂、二「•.母二-.•- : - ... 丁 5）
①【,上_.7这是一个k 行k-1列的Hessenberg 矩阵，如果将 Avk 用Aiv （i>=0）表示出来，再代
hi-lJc-2 1=1
-
h'1,1 h'1,2 III h\,k v Av |
II A k
」v]=Iv 1 V 2 川 V k ]
0 +
+ h'2,2 III + b + B h'2,k
1- ►
1
0 0 III h'k,k
②其实该方法每一步都有用到 Gram-Schmidt 正交化方法；但与 Schmidt 正交化不同的地方 是，该方法用［替换向量组（H ）中的 「得：,■,［，駅毀•，新得到的向量组与前面的 向量组都是等价的，如此反复，便得到一组标准正交基
Schmidt 正交化方法 中
没有本题算法的替代过程。

2.3实习题目3 题目
已知A€Rm 为对称正定矩阵，二次函数
林X ）
1
i
i n n
x x T
Ax-b T x
Ax,x I ib,x :
2 2
2 y jd ：
在x 处，若a€R ,设u = （U ：, u 2, ...,u 门）丁€卍为一方向向量， （1 ）请你计算如下的导数：
——;：
X 上：心 （3-
2 ）
d «
并利用内积表示该导数。

（2）当u 为单位向量时，请你根据 Cauchy - Schwarz 不等式，计算以上导数绝对值的
最大值。

求解过程: 1)
对——切 x,u • R n , ：£ e
R
(x 亠很⑴ (A(x 亠MU ), x 亠 u) - (b, x 亠 u)
2
1 : 2
(Ax, x) —(b,x)亠：：£(Ax, u) - ： (b, u) (Au, u) 2 2
••再求一阶导，得 。

丘卜討））-加 • if （做吊）=（相艾"疵 心我）
eta
2）根据 Cauchy-Schwarz 不等式 川-卜限
Schmidt 正交化方法产生的上三角矩阵
n
ajXXj-' bx
i 二
(3-1 )
=(x) ： (Ax —b,u) (Au, u)
我们有(訂2亠越)—MiQ — 亠…; —.；• 卜，
又因为•是单位向量，「工「代入得 •••( . 一 1 •. 一 ：’ .)： ./(卜农,:计 由此可知，该函数在某一固定点沿梯度方向变化最快。

证毕 2.4实习题
目4
题目 已知A € R n Xn
为对称正定矩阵，定义
1
X A 二 Ax,x°
设A 的特征值 入排序为0 w 入w ... ＜乩p(t)是t 的一个多项式, h(A )X ^max
P " H X |A ，P x W R
n
(4-2
)
成立。

(提示：利用对称正定线性方程组的特征向量组与标准正交基之间的关系进行证明) 求解过程： 设 I 、：：：：；.‘：‘.■；「 是 A 的 对 应 于
、
「丁订苛匸—扌山=口卜「卡'I-'- ：■-,有菸
V 立,因为 pU)y ： = pUt)yr ，进
而有
n TI
(pU)xp(4)x)4=片⑷如(如(0潭(人)娜①丽⑷旳
(=1
i=i
戈做弊，鴻除曲沁①松曲品驚庞佻=.一一馅心厂⑴|习.；：.
(
4-1
)
请你证明如下的不等式：
即:泸肩协就桃圧咏轴J 関站闊b
2.5实习题目5 题目
已知Ax = b 为对称正定线性方程组，设 p, q € R n
,我们说p 和q 是A 共轭的，是指它
们满足（p, q ）A = p T
Aq = 0。

设p i ,…,p n-k 是R n
中n - k 个线性无关向量， 讯…,沁是给定的n - k 个实数。

我们称满
足
P T
X Y , i =1,2,11）,n-k
（
5-1
）
的X 的全体所构成的 R n
中的点集为一个k 维超平面，记作
■ :k : p T
x =订，i =1,2，川，n-k
（
5-2
）
其中的P i 称作n 的法向量。

设n 是任一以P I , ..., p n-k 为法向的k 维超平面，则称由 A -1
p 1, ..., A -1
p n-k 所张成的超平面
玄上：X •川 VnJ^A’P n 丄
（
5-3
）
为■:k 的共轭超平面，其中 X * = A -1
b 是已知方程组的解。

显然， ？上的维数为n - k 。

请你证明：
二k 内的向量与 W 上内的向量是互相共轭的， 而且若二k 是由U 1, ..., u k 张成的，则Au 1, ..., Au k 就是2上的法向量。

（提示：将线性组合利用矩阵进行表示）
求解过程：
任取？n
_k 中向量q ，则q = q 1 - q 2, q 1, q 2《，于是有:
、厂 n 」
P i - |x *+2： Ri （2）A 」P i
丿I 7
n -k
def n -k
def
=Z （刖）-郎2）沪和冃人和=A 〜p0
i =1
i
A
任取n 内的向量r ，由题意可知：
P T
r =0, i =1,2，川,n -k
记p =(P 1, p 2, ..., p n-k )'则p T
r=0n-k,1，再根据本题内积定义，有：
r, q 二 r T
Aq 二 r T A A'-
_r T p- - p T r : =0
所以二k 内的向量与 内的向量是互相共轭的。

进而可以推导出：
n -k
q=q —q 2 =|x * 吃 g A 」
I 7
def n -k
(5-4)
(5-5)
(5-6)
(5-7)
Am 爲=Ui T A T q Aq 二⑴,q ]=。

，i =1,2」|| ,k 得证AU l,…，Au k就是：?n上的法向量。

2.6实习题目6
题目
请你编写一个程序，对一个长方体进行矩形块剖分，按照先x方向，然后y方向、最后z方向的顺序对这个剖分结构进行编号，并能提取任一平行坐标平面的截面及所对应的矩形块编号。

求解过程或算法：
function [A,g] = Split (x,y,z,w,q) %A是矩阵g是剖分矩阵，w,q是考察魔方w方向的第q层的标号,g是剖分的,
A是二维表示三维的标号
A=on es(x,y*z);
for i=1:x
for j=1:y*z
A(i,j)=1+(i-1)+x*(j-1);
end
end
D=A;
switch w
case 1 %对y进行剖面
g=[];
u=[];
for j=1:z
n=j;
u=D(q,(n-1)*y+1:n*y);
u=u;
g=[g;u];
end
case 2 %对x进行剖面
g=[];
u=[];
for n=1:z
u=D(:,(n-1)*y+q);
g=[g,u];
end
case 3 %对z进行剖面
g=D(1:x,(q-1)*y+1:q*y);
end
end
程序结果：
>> [A,g]= Split (3, 3, 3, 3, 1)
1
4 7 10 13 16 19 22 2
5 2 5 8 11 14 17 20 23 2
6 3
6
9
12 15 18
21 24
27
g =
1 4 7 2
5
8
3 69 2.7实习题目7 题目
已知无补给的承压水完整井非稳定流
Theis 井流公式：
Q f
2
缶-0.577216-1nu+u-J
4 兀T i
4
其中:
2
r s u 二 4Tt
Q —抽水井的流量；r —观测井与抽水井的距离； t —抽水开始到计算时刻的时间（抽水
时间）；S cal —水位降深的计算值。

已知
S obs 为水位降深的实测值，如下表：
各时间、井距的水位降升数据
(7-1)
(7-2)
初始值取：Q o = 0.0167m3/s, T o = 0.001m2/s, s o = 0.0001。

请你根据表上的数据计算出Theis井流公式的参数s （贮水系数）和 T （导水系数）。

求解过程或算法：
本题可采用最小二乘法的思想，寻找使得残差平方和最小的S1：赋初值S2:给定最大迭代次数和误差限S4:计算残差和S5 :计算梯度
程序：
function [zO,varParameter] = cal_coefficient(g)
global A varRadius varTime varDropDepth varQuantity
%设置全局变量，注意顺序
A=xlsread( 'pump' ); %读取并导入excel的数据varRadius=A(:,3);
varTime=A(:,1);
varDropDepth=A(:,2);
varQuantity=0.0167;
[m,n]=size(varRadius');
zSingle=[];
varParameter=[0.001,0.0001]; % 初值
Fri=zeros(1,2);
upper=0;
lower=0;
z0=Residual(n,varParameter); T、s。

具体算法如下：
S3：计算变化值
S6：计算二阶导数
S7:计算海塞矩阵S8:计算修正值
S10：判断结束迭代的条件若满足，则结束S3继续迭代。

S9:计算新目标值
T i和s i即为所求。

否则 T = T 1，s = s 1，转
for w=1:g
for varIndex=1:2
dltPara=0.001*varParameter(varIndex);
if varParameter(varIndex)==0
dltPara=0.0001
end
varParameter(varlndex)=varParameter(varlndex)+dltPara; % 力口上一个A0
zSingle(varlndex)=Residual(n,varParameter); %调用ObjFun 计算Z(x,(+ △》
varParameter(varlndex)=varParameter(varlndex)-dltPara; %还原一个参数值，循环计算另一个end for varSequence = 1:2
dltParaseq=0.001*varParameter(varSequence);
if varParameter(varSequence)==0
dltParaseq=0.0001;
end
varParameter(varSequence) = varParameter(varSequence) + dltParaseq;
for varlndex = 1:2
dltParaind = 0.001 * varParameter(varlndex);
if varParameter(varlndex) ==0
dltParaind=0.0001;
end
varParameter(varlndex)=varParameter(varlndex) + dltParaind;
zMultipl e(varSequence,varlndex)=Residual(n,varParameter); %计算海赛矩阵
varParameter(varlndex)=varParameter(varlndex) - dltParaind;
end
varParameter(varSequence)=varParameter(varSequence) - dltParaseq; end% 计算二阶导数矩阵的前期准
备工作
zMultipl e
for varlndex =1:2
Fri(varlndex)=(zSingl e(varlndex)-z0)/(0.001*varParameter(varlndex)); %计算一阶导数
for Seq =1:2 %计算二阶导数矩阵
Sec(varlndex, Seq) = (zMultipl e(varlndex, Seq)-zSingl e(varlndex)-zSingl e(Seq) +z0)/(0.001 *
0.001*varParameter(varlndex)*varParameter(Seq));
end
end
% 计算阿拉法 a 的值
for varlndex = 1:2
g=Fri(varlndex)*Fri(varlndex);
upper=upper+g;
for Seq =1:2
lower=l ower+Fri(varlndex)*Sec(varlndex, Seq)*Fri(Seq); end
end
alf=upper/l ower;
%计算参数的修正值
for varIndex=1:2
varParameter(varlndex)=varParameter(varlndex)-alf*Fri(varlndex);
end
%用新的参数算出新的z值
zNew = Residual(n,varParameter);
%跳岀循环条件
if (abs(z0-zNew)/z0<0.000001) % 设定相对误差限
break;
else
z0 = zNew;
end
end
end function [Sum] =Residual(n,w)
global A varRadius varTime varDropDepth varQuantity %全局变量
varParameter=w;Sum=0; %待求参数s,T组成的二维向量for varlndex =1:n % 向量下标
Uvalue = varRadius(varlndex)*varRadius(varlndex)*varParameter(2)/(4 * varParameter(l) * varTime(varlndex));
Sqt = varQuantity /(4*pi*varParameter(1))*(-0.577216-l og(Uvalue) +
Uvalue-Uvalue*Uvalue/4)-varDropDepth(varlndex);
Sum = Sum+Sqt*Sqt; % 误差求和
end
end
程序结果：
>> cal_coefficie nt(100000)
得到返回参数导水系数T=0.0029二一：，贮水系数s=0.000105，minZ=3.0893 。

2.8实习题目8
题目
在简单蒸馏釜内蒸馏 1000kg含C2H5OH（乙醇）质量分数为60%和出0质量分数为40% 组成的混合液。

蒸馏结束时，残液中含C2H5OH的质量分数为5%,请你利用简单蒸馏的雷
利公式
,F X F dx ln
W *w y —x(8-1)
求残液的质量。

其中，F 为原料液量，W 为残液量，X F 为原料液组成，X W 为残液组成。

该体系的气液平衡数据如表 1所示，其中x 为液相中C 2H 5OH 的质量分数，y 为气相中C 2H 5OH 的质量分数。

求解过程或算法:
根据题意分析可以知道 F=1000kg ，轧：址云一工迂，又因为y 和x 存在某种函数 关系，所以我们不能直接积分。

又发现求积节点等距分布， 所以我考虑利用数值分析中的复
化求积公式。

将1/( y - x )记为f (x )，雷利公式右端的积分项记为 I ，对前11个数据点利
用复化Simpson 公式对最后两个数据点利用梯形公式， 求得数值积分值。

其中步长h 取0.05。

然后用如下表达式计算残液量
W
(8-2)
复化梯形公式:
复化Simpson 公式(此时 h = 0.1)
程序：
function run x=[0.05:0.05:0.6];
y=[3.22 2.4 2.22 2.2 2.27 2.44 2.64 2.9 3.29 3.74 4.38 5.29]; [m,n]=size(y); k=1; q=k+2; sum=0; for i=1:n
h=(1/6*y(k)+4/6*y(k+1)+1/6*y(q))*(0.15-0.05); sum=sum+h; %对前五个区间段利用辛普森公式求和 k=k+2;
ln F
_L 二 e
h X W
X F
2
n 4
亠二 hf x
i
=1
(8-4)
LUX W 丄 2h^
6
3 i ^0
f
x
i 1/2
(8-5)
q=k+2; if(q>n) break; end end
sum=sum+(0.6-0.55)*1/2*(y(n-1)+y(n));
% 利用梯形公式处理
q=exp(l og(IOOO)-sum) %对原积分公式进行变形求解
W
程序结果：
>> run
q = 196.3709 >>
结论：算出残液量为 196.3709kg 。

2.9实习题目9 题目
某连续反应过程
k 2
k 3
k 4
A > B_ C_ D_ E
体系中各物质的动力学方程分别为
已知 k = 0.01s -1
, k 2 = 0.20s -1
，k 3 = 0.10s -1，k 4 = 0.05s -1。

反应开始时只有
A 物质存在,
其浓度为1mol/L ，试求前25s 中每隔5s 时各物质的浓度(取计算步长 0.5s)。

求解过程或算法：
本题是一阶线性常微分方程组的数值求解问题，利用四阶 R-K 公式的向量形式确定具
体的迭代公式，并表成简洁的矩阵形式， 利用MATLAB 进行编程计算，计算结果为矩阵 Y ，
Y 的第i 行为第5i s 时的各物质的浓度。

初值取C °=(1,0, 0, 0, 0)。

由四阶R-K 方法，可得迭代格式：
(9-1)
dc A dt dC B dt =k^A - k ?C B
dC C dt
-k 2C
^ _ k 3C
C
(9-2)
罟二 k3CC
芽
=
k4CD
-k 4C
D
程序：
function [C] = RKmatrix(T) a1=0.01; a2=0.20; a3=0.10;
a4=0.05; %参数值 h=0.5; % 步长 0.5 n=T/h+1; A=[-a1 0 0 0 0
a1 -a2 0 0 0 0 a2 -a3 0 0 0 0 a3-a4 0 0 0 0 a4 0];
C=zeros(5,n); C(1,1)=1; for i=1:n-1
k1=h*A*C(:,i); k2=h*A*(C(:,i)+k1/2); k3=h*A*(C(:,i)+k2/2); k4=h*A*(C(:,i)+k3);
C(:,i+1)=C(:,i)+(k1+2*k2+2*k3+k4)/6; %四阶龙格库塔方法的矩阵形式
if(rem(i,10)==0)
C(:,i+1) %判断余数，并输出每间隔
5s 的数据
end end
程序结果：
输入>> RKmatrix(25)得到各物质浓度
前25s 中每隔5s 时各物质的浓度如下表所示。

各时刻各物质浓度
2.10实习题目10
其中, C ni 二5 £ g l 2g 2 2g 3 g 4
(9-4qizh
g i 二 KC n , g 2 二 K ! C n -2 ,
g^K C n 号，
g^K C n hg 3
(9-5)
题目
对套管-地层系统在均匀地应力条件下的套管载荷进行弹塑性理论分析时，采用 Mohr-Coulomb 准则，给出了受内外压厚壁筒弹、塑性区的位移和应力分布的解析解，由该 解析解可直接退化得到 Tresca 准则条件下厚壁筒的位移和应力分布。

将套管
-地层系统看作
是岩石和钢材两种不同材料厚壁筒的组合，在解析解的基础上，可以推导出套管 -地层系统
仅地层进入塑性和仅套管进入塑性两种情况下的套管载荷和极限地应力表达式。

将地层简化为弹性区域，套管视为弹塑性区域（ Tresca 准则），考虑套管进入塑性状态
下的套管载荷和地应力的关系，得出地应力
b 为
其中r p 为塑性半径（m ）; &为套管载荷（MPa ）; a °, a 分别为套管内、外径（m ）; E s , E c 分别为地层和套管的弹性模量 （MPa ） ; V s , V c 分别为地层和套管的泊松比 （无量纲）; b 为套管外壁应力（MPa ）。

问题：若已知系统的参数 a 。

= 62.13mm , a = 69.85mm , E s = 40GPa, E c = 210GPa,V s = 0.18 , V c = 0.3 , b = 770MPa ，为了保持系统状态，取 b [107.7, 128.57]，当 b = 120MPa 时，请计 算套管的塑性半径 ％。

求解过程或算法：
本题就是非线性方程的根的数值求解问题, 易地求出套管的塑形半径 r p 程序：
function [x,k]=cal culate(x,N,e) k=1;a0=62.13; a=69.85;Es=40; Ec=210;vs=0.18;
vc=0.3;b=770; % 初始条件 C1= b/(2*(1-VS));
C2=((1+VC )*E S )/((1+VS )*E C ); C3=(1-2*VC );
C4=(1-2*vc)/2; whil e k<=N
x0=x;
f=120-(C1*l og(x/a0)+1/2*C1*(1-(x*x)/(a*a))+C1*C2*C3*l %f(x)=0
df=-(C1/x-C1*x/(a*a)+C1*C2*C3/x+C1*C2*x/(a*a)); x=x-f/df;
%牛顿迭代法
if abs(x-x0)<e % 终止条件
break
end k=k+1; end
程序结果：
由题目已知参数取初值为 66,最大迭代次数1000,误差限0.00001.
>> [x,k]=calculate(66,1000,0.000001)
•丄二—V c ln 电上空2
1 V s E c
c
a o 2 2 a 2
(10-1)
利用经典的牛顿迭代法进行编程可以较为容
og(x/a0)+C1*C2*C4+(C1*C2)*(x*x)/(2*a*a));
%f 的导数
2 r
p
a
66.1477
>>
用牛顿迭代法求得套管的塑性半径r p = 66.1477mm，迭代次数为3次。

2.11实习题目11
题目：
在半径为R = 2m的圆筒形贮槽中，开始时加水至M = 5m，然后用半径门=2 X10-2m的
给水管以稳定的流速V1 = 0.7m/s的速度向槽内加水，同时，由位于槽底部半径为「2 = 7.6X10-2m的排水管排水。

如果不考虑排水管的压头损失，试求开始排水后前4分钟内每分
钟贮槽水位的高度y (m)。

求解过程或算法：
其数值解可用4阶R-K方法计算，取步长 h = 1s,则迭代格式如下：
y n 1 二y n 一K 2k2 2k3 k4 (11-2
)
6
其中，
k^ f y n , k^ f y n hk / 2 , k3 = f y n hk2/2 , k^ f y n hk? (11-3)
程序：
function [ y ] = RK( n )
R=2;r1=2*0.01;r2=7.6*0.01;
v1=0.7;M=5;g=9.8;h=1; % 一些初始条件
y=zeros(1,60*n+1);
y(1)=5; %初边值
t=0:h:60*n;
for n=2:h:60*n+1
k1= (0.7*r1*r1-r2*r2*sqrt(y(n-1)*2*g))/(R*R);
k2=(0.7*r1*r1-r2*r2*sqrt((y(n-1)+1/2*h*k1)*2*g))/(R*R);
k3=(0.7*r1*r1-r2*r2*sqrt((y(n-1)+1/2*h*k2)*2*g))/(R*R);
k4=(0.7*r1*r1-r2*r2*sqrt((y(n-1)+h*k3)*2*g))/(R*R);
y(n)=y(n-1)+h/6*(k1+2*k2+2*k3+k4);
end plot(t,y,'r-')
title('桶内水位变化'); xlabel('时间 t'); ylabel('高度 y'); end
程序结果：
在 comma nd window 里输入 >> RK(4)
得到结果，开始排水后前4分钟内每分钟贮槽内水位的高度分别是 4.1831m 、 2.7688m 、2.1712m 。

过程如下图所示。

Figure 1 — □ X
文忤(£)駅捐宜看(旳插入皿工具
(D 宴面血)闔口的帮助(H 】
Odd 2 4 7、彰越哽山•思口因■ O
%四阶龙格库塔方法
3.4394m 、。