中国石油大学计算实习报告.docx
- 文档编号:12271809
- 上传时间:2023-06-05
- 格式:DOCX
- 页数:32
- 大小:216.77KB
中国石油大学计算实习报告.docx
《中国石油大学计算实习报告.docx》由会员分享,可在线阅读,更多相关《中国石油大学计算实习报告.docx(32页珍藏版)》请在冰点文库上搜索。
中国石油大学计算实习报告
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
0:
:
x:
:
1,t0
汀一:
:
T
一;x2
求解过程或算法:
首先参考《数值分析》P67的方法,编写了列主元高斯消去法的matlab程序:
function[X]=Gauss(B)
[u,v]=size(B);%求出输入矩阵B的维度
m=u;
X=zeros(m,1)%类似一个初始化,建立一个m行1列的0矩阵
forn=1:
1:
m-1%一个循环嵌套,步长为1,从1到m-1行(因为最后一行m行对角线上元素不取)
u=B(n:
m,n);%建立一个矩阵n行到m行,位置第n列的一个(m-n)行1列的矩阵[c,i]=max(u);%求里面的最大元素的位置
ifc==0%如果列主元为0,贝UdetA=0
error('error');
else
ifi+n-1~=n%如果不在对角线上
t=B(n,:
);%取一个中间变量,最大元所在n行的一整行元素
B(n,:
)=B(i+n-1,:
);
B(i+n-1,:
)=t;%把列主元所在的第n行和第i+n-1行整行对换
end
end
forj=n+1:
1:
m
B(j,:
)=B(j,:
)-B(j,n)./B(n,n)*B(n,:
);%对角元以下消元计算,变成一个上三角矩阵
endend
A=B(:
1:
m);%1到m列的数据
d=B(:
m+1);%增广矩阵最后一列,即b
fork=m:
-1:
1%倒序
z=0;
forl=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=length(x);
z=zeros(T+1,m);%温度矩阵
u=sin(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函数)
forn=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
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');
程序结果:
在commandwindow里运行CN(1,20),我们可以得到如下两张结果图:
图1图2
结论:
由图1可以看出,在网格比/二:
[时,得出的数值解与精确解非常近似,这是因为C-N
格式是无条件稳定的。
2.2实习题目2
题目设v是非零n维向量,A是nxn可逆矩阵,与A相关的向量组为
v,Av,A2v,III,Ak4v(2-1)
线性无关,记Kk=span{v,Av,A2v,...,Ak-1v}。
。
以k=4为例,如果按如下的计算过程:
1将v单位化:
令V1=v/||v||,然后用v1替换原向量组中的v,得
23
v1,Avi,Av-!
Av1(2-2)
该向量组显然与原向量组等价。
2利用Gram-Schmidt正交化方法,从Av1中减去它在v1上的投影向量,得
(2-3)
v2=Aw-Av1,viv
单位化后得V2=V2/||V2|。
用V替换向量组(2-2)中的AV1,得:
Vi,V2,Av,,A2V2
(1)请你证明:
向量组(I)和(II)等价,即它们可以相互线性表示。
(2)请你按照①、②的逻辑过程,继续对向量组进行执行计算、替换操作,写出计算过程,最终得到如下的规范化正交向量组(标准正交基)
(2-4)
Vi,V2,V3,V4
hi,j二AVj,Vi,hji,j=Vji
(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正交化方法,从空%:
中减去它在馬「:
底上的投影向量,得
.—单位化后得&:
;-:
川〔;|。
用•一替换向量组(I)中
的血雄,得:
一.殴飞麹.
利用Gram-Schmidt正交化方法,从」i中减去它在说|:
:
险險上的投影向量,得
二一.二]:
.:
单位化后得用_替换向
量组(I)中的姒;,得:
3)算法:
1.开始:
._=』|忌||
2•迭代
forj=1k-1do
hjj=巧=1,2j
i
町+1=阿_2如円
i
End4)将第二问的结果写成矩阵形式,结合第三问的结果:
■
观察算法中式子的'以及占曲隔;.得到
矽为上下两块©)心釉吩别梯怖和最后-彳亍
叫7陋,Hk为k-1阶方阵,且胪的前k-1个元素都为零
^k-i>fc-v
•••归纳出的以上算法的矩阵形式为厂、二「•.母二-.•-:
-...丁5)①【,上_.7这是一个k行k-1列的Hessenberg矩阵,如果将Avk用Aiv(i>=0)表示出来,再代
-
h'1,1
h'1,2III
h\,k
vAv|
IIAk」v]=Iv1V2川Vk]
0
+
+
h'2,2III
+b
+B
h'2,k
1-
►
1
0
0III
h'k,k
Schmidt正交化方法产生的上三角矩阵
②其实该方法每一步都有用到Gram-Schmidt正交化方法;但与Schmidt正交化不同的地方是,该方法用[替换向量组(H)中的「得:
■,[,駅毀•,新得到的向量组与前面的
向量组都是等价的,如此反复,便得到一组标准正交基Schmidt正交化方法中
没有本题算法的替代过程。
2.3实习题目3
题目
已知A€Rm为对称正定矩阵,二次函数林X)
1
iinnxxTAx-bTxAx,xIib,x:
222yjd:
在x处,若a€R,设u=(U:
u2,...,u门)丁€卍为一方向向量,
(1)请你计算如下的导数:
——;:
X上:
心(3-2)
d«
并利用内积表示该导数。
(2)当u为单位向量时,请你根据Cauchy-Schwarz不等式,计算以上导数绝对值的
最大值。
求解过程:
1)
对——切x,u•Rn,:
£eR
(x亠很⑴(A(x亠mu),x亠u)-(b,x亠u)
2
1:
2
(Ax,x)—(b,x)亠:
:
£(Ax,u)-:
(b,u)(Au,u)
22
••再求一阶导,得。
丘卜討))-加•if(做吊)=(相艾"疵心我)
eta
2)根据Cauchy-Schwarz不等式川-卜限
我们有(訂2亠越)—MiQ—亠„;—.;•卜,
又因为•是单位向量,「工「代入得•••(.一1•.一:
’.):
./(卜农,:
计
由此可知,该函数在某一固定点沿梯度方向变化最快。
证毕2.4实习题目4
题目已知A€RnXn为对称正定矩阵,定义
1
(4-1)
请你证明如下的不等式:
Xa二Ax,x°
设A的特征值入排序为0w入w...<乩p(t)是t的一个多项式,
h(A)X^maxP"Hx|a,PxWRn(4-2)
成立。
(提示:
利用对称正定线性方程组的特征向量组与标准正交基之间的关系进行证明)求解过程:
设I、:
:
:
:
;.‘:
‘.■;「是A的对应于、
「丁订苛匸—扌山=口卜「卡'I-'-:
■-,有菸V立,因为pU)y:
=pUt)yr,进而有
nti
(pU)xp(4)x)4=片⑷如(如(0潭(人)娜①丽⑷旳
(=1i=i
戈做弊,鴻除曲沁①松曲品驚庞佻=.一一馅心厂⑴|习.;:
.
即:
泸肩协就桃圧咏轴J関站闊b
2.5实习题目5
题目
已知Ax=b为对称正定线性方程组,设p,q€Rn,我们说p和q是A共轭的,是指它
们满足(p,q)A=pTAq=0。
设pi,…,pn-k是Rn中n-k个线性无关向量,讯…,沁是给定的n-k个实数。
我们称满
足
PtxY,i=1,2,11),n-k(5-1)
的X的全体所构成的Rn中的点集为一个k维超平面,记作
■:
k:
pTx=订,i=1,2,川,n-k(5-2)
其中的Pi称作n的法向量。
设n是任一以Pi,...,pn-k为法向的k维超平面,则称由A-1p1,...,A-1pn-k所张成的超平面
玄上:
X•川VnJ^A’Pn丄(5-3)
为■:
k的共轭超平面,其中x*=A-1b是已知方程组的解。
显然,?
上的维数为n-k。
请你证明:
二k内的向量与W上内的向量是互相共轭的,而且若二k是由U1,...,uk张成的,则Au1,...,
Auk就是2上的法向量。
(提示:
将线性组合利用矩阵进行表示)
求解过程:
任取?
n_k中向量q,则q=q1-q2,q1,q2《,于是有:
、厂n」
Pi-|x*+2:
Ri
(2)A」Pi
(5-5)
(5-6)
(5-7)
丿I7
n-kdefn-kdef
=Z(刖)-郎2)沪和冃人和=A〜p0
i=1iA
任取n内的向量r,由题意可知:
PTr=0,i=1,2,川,n-k
记p=(P1,p2,...,pn-k)'则pTr=0n-k,1,再根据本题内积定义,有:
r,q二rTAq二rTAA'-_rTp--pTr:
=0
所以二k内的向量与内的向量是互相共轭的。
进而可以推导出:
Am爲=UiTATqAq二⑴,q]=。
,i=1,2」||,k
得证AUl,…,Auk就是:
?
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=ones(x,y*z);
fori=1:
x
forj=1:
y*z
A(i,j)=1+(i-1)+x*(j-1);
end
end
D=A;
switchw
case1%对y进行剖面
g=[];
u=[];
forj=1:
z
n=j;
u=D(q,(n-1)*y+1:
n*y);
u=u;
g=[g;u];
end
case2%对x进行剖面
g=[];
u=[];
forn=1:
z
u=D(:
(n-1)*y+q);
g=[g,u];
end
case3%对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
25
2
5
8
11
14
17
20
23
26
3
6
9
12
15
18
21
24
27
g=
147
258
3692.7实习题目7
题目
已知无补给的承压水完整井非稳定流Theis井流公式:
Qf2
(7-1)
(7-2)
缶-0.577216-1nu+u-J
4兀Ti4
其中:
2rsu二
4Tt
Q—抽水井的流量;r—观测井与抽水井的距离;t—抽水开始到计算时刻的时间(抽水
时间);Scal—水位降深的计算值。
已知Sobs为水位降深的实测值,如下表:
各时间、井距的水位降升数据
井号
抽水时间(s)
降深(m)
井距(m)
井号
抽水时间(s)
降深(m)
井距(m)
1
600
0.73
43
2
600
0.16
140
1
1200
1.28
43
2
1200
0.48
140
1
1800
1.53
43
2
1800
0.54
140
1
2400
1.72
43
2
2400
0.65
140
1
3600
1.96
43
2
3600
0.75
140
1
4800
2.14
43
2
4800
1
140
1
6000
2.28
43
2
6000
1.12
140
1
7200
2.39
43
2
7200
1.22
140
1
9000
2.54
43
2
9000
1.36
140
1
12600
2.77
43
2
12600
1.55
140
1
16200
2.99
43
2
16200
1.7
140
1
19800
3.1
43
2
19800
1.83
140
1
24000
3.2
43
2
24000
1.89
140
1
27000
3.26
43
2
27000
1.98
140
1
38700
3.47
43
2
38700
2.17
140
1
52200
3.68
43
2
52200
2.38
140
1
59400
3.77
43
2
59400
2.46
140
1
71100
3.85
43
2
71100
2.54
140
井号
抽水时间(s)
降深(m)
井距(m)
井号
抽水时间(s)
降深(m)
井距(m)
3
600
0.04
510.00
3
59400
1.24
510.00
3
2400
0.06
510.00
3
71100
1.35
510.00
3
3600
0.20
510.00
4
4800
0.04
780.00
3
4800
0.20
510.00
4
7200
0.08
780.00
3
6000
0.20
510.00
4
9000
0.09
780.00
3
7200
0.21
510.00
4
12600
0.16
780.00
3
9000
0.24
510.00
4
16200
0.25
780.00
3
12600
0.40
510.00
4
19800
0.34
780.00
3
16200
0.53
510.00
4
24000
0.42
780.00
3
19800
0.63
510.00
4
27000
0.50
780.00
3
24000
0.65
510.00
4
38700
0.71
780.00
3
27000
0.73
510.00
4
52200
0.87
780.00
3
38700
0.93
510.00
4
59400
0.96
780.00
3
52200
1.14
510.00
4
71100
1.06
780.00
初始值取:
Qo=0.0167m3/s,To=0.001m2/s,so=0.0001。
请你根据表上的数据计算出
Theis井流公式的参数s(贮水系数)和T(导水系数)。
求解过程或算法:
本题可采用最小二乘法的思想,寻找使得残差平方和最小的
S1:
赋初值S2:
给定最大迭代次数和误差限
S4:
计算残差和S5:
计算梯度
S7:
计算海塞矩阵S8:
计算修正值
S10:
判断结束迭代的条件若满足,则结束
S3继续迭代。
程序:
function[zO,varParameter]=cal_coefficient(g)
globalAvarRadiusvarTimevarDropDepthvarQuantity
%设置全局变量,注意顺序
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);
forw=1:
g
forvarIndex=1:
2
dltPara=0.001*varParameter(varIndex);
ifvarParameter(varIndex)==0
dltPara=0.0001
end
varParameter(varlndex)=varParameter(varlndex)+dltPara;%力口上一个A0
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 中国 石油大学 计算 实习 报告