2022年西安交通大学2021年计算方法A上机大作业 .pdf

《2022年西安交通大学2021年计算方法A上机大作业 .pdf》由会员分享，可在线阅读，更多相关《2022年西安交通大学2021年计算方法A上机大作业 .pdf（13页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、计算方法 A上机大作业张晓璐硕 4011班学号： 3114009097 1. 共轭梯度法求解线性方程组算法原理： 由定理 3.4.1 可知系数矩阵 A是对称正定矩阵的线性方程组Ax=b的解与求解二次函数1( )2TTf xx Axb x极小点具有等价性，所以可以利用共轭梯度法求解1( )2TTf xx Axb x的极小点来达到求解Ax=b的目的。共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式：(1)()( )kkkkxxd产生的迭代序列(1)(2)(3)xxx，.在无舍入误差假定下，最多经过n 次迭代，就可求得( )f x的最小值，也就是方程Ax=b的解。首先导出最佳步长

2、k的计算式。假设迭代点()kx和搜索方向()kd已经给定，便可以通过( )( )( )()kkf xd的极小化( )()min( )()kkf xd来求得，根据多元复合函数的求导法则得: ()()( )( )()kkTkf xdd令( )0，得到 : ( )( )( )( )k Tkkk TkrddAd，其中( )( )kkrbAx然后确定搜索方向( )kd。给定初始向量(0)x后，由于负梯度方向是函数下降最快的方向，故第一次迭代取搜索方向(0)(0)(0)(0)()drf xbAx。令(1)(0)00 xxd其中(0)(0)0(0)(0)TTrddAd。第二次迭代时，从(1)x出发的搜索方向

3、不再取(1)r，而是选取(1)(1)(0)0drd，使得(1)d与(0)d是关于矩阵 A的共轭向量， 由此可求得参数0: 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 13 页(1)(0)0(0)(0)TTrAddAd然后从(1)x出发，沿(1)d进行搜索得到(2)(1)(1)1xxd设已经求出(1)( )( )kkkkxxd，计算(1)(1)kkrbAx。令(1)(1)( )kkkkdrd，选取k，使得(1)kd和( )kd是关于 A 的共轭向量，可得: (1)()( )( )kTkkk TkrAddAd具体编程计算过程如下：（1）

4、 给定初始近似向量(0)x以及精度0；（2） 计算(0)(0)rbAx，取(0)(0)dr；（3） For k=0 to n-1 do （i ）( )()0( )( )k Tkk TkrddAd；（ii ）(1)()( )kkkkxxd；（iii）(1)(1)kkrbAx；（iv ）若(1)kr或k=n-1，则输出近似解(1)kx，停止；否则，转（ v） ；（v）2(1)22( )2kkkrr；（vi ）(1)(1)()kkkkdrd；End do 程序框图：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 13 页开始输入 , , ,

5、(0)x(0)(0)rbAx(0)(0)dr( )( )( )( )k Tkkk TkrddAd(1)( )()kkkkxxd(1)(1)kkrbAx0k1kkAb( ,1)nsize A或(1)kr1knT输出近似解(1)kxF结束2(1)22()2kkkrr(1)(1)( )kkkkdrd程 序 使 用 说 明 ： 本 共 轭 梯 度 法 求 解 线 性 方 程 的 程 序 是 一 个 函 数Gongetidu2(A,b,x0,epsa)，在求解线性方程组Ax=b（A 是对称正定矩阵）的时候 ， 首 先 给 定 初 始 向 量x0和 误 差epsa ， 然 后 直 接 调 用 该 函 数G

6、ongetidu2(A,b,x0,epsa)即可，其中 A,b,x0 和 epsa 都是可变的，虽然该函数是通用的，但是对于矩阵 A和向量 b 的输入，需要使用者根据 A和 b 的特点自行输入。算例 3.2 计算结果：对题 113 页 3.2 ，首先编程将矩阵A,b,x0,epsa读入系统，然后再调用函数 Gongetidu2(A,b,x0,epsa)；程序如下：clear A b x0 %程序运行前首先清除A，b 和 X0的数值，以免影响计算clc n=input(请输入对称正定矩阵A的阶数 n=); epsa=input(请输入误差 =); for k=1:(n-1) %读取题目 3.2

7、的矩阵 A A(k,k)=-2; A(k,k+1)=1; A(k+1,k)=1; end 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 13 页A(n,n)=-2; A; b(1,1)=-1; %读取题目 3.2 的矩阵 b b(n,1)=-1; b; x0(1,1)=input(假定初始向量各元素相同，均为：); %给题目 3.2 的迭代过程赋初值for kk=2:n x0(kk,1)=x0(1,1); end x=Gongetidu2(A,b,x0,epsa) %调用共轭梯度法求线性方程组的函数function x=Gongeti

8、du2(A,b,x0,epsa) n=size(A,1); x=x0; r=b-A*x; d=r; for k=0:(n-1) alpha=(r*r)/(d*A*d); x=x+alpha*d; r2=b-A*x; if (norm(r2)=epsa)|(k=n-1) x; break; end beta=norm(r2)2/norm(r)2; d=r2+beta*d; r=r2; end 计算结果：当 n=100时，x=1;1;1;1;1;1;1;.1;1;1;1 当 n=200时，x=1;1;1;1;1;1;1;.1;1;1;1;1;1;1 当 n=400时，x=1;1;1;1;1;1;1

9、;.1;1;1;1;1;1;1;1;1;1 2. 三次样条差值算法原理（三次样条插值函数的导出） ：(i).导出在子区间1,iixx上的 S(x) 的表达式由于 S(x) 的二阶导数连续，设 S(x) 再节点ix处的二阶导数值 S(xi)=Mi ，其中 Mi 为未知的待定参数。由S(x) 是分段的三次多项式知，S(x) 是分段线性函数， S(x) 在子区间1,iixx上可表示为精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 13 页1111111( ),iiiiiiiiiiiiiiiixxxxSxMMxxxxxxxxMMxxxhh其中

10、hi=xi-x(i-1),对上式两次积分得到331111( )66()(),iiiiiiiiiiiixxxxS xMMhhb xxcxxxxx由插值条件11(),()iiiiS xyS xy得到2211() /,() /66iiiiiiiiiihhbyMh cyMh将iibc和代入( )S x可得3321111211( )()/()666()/(),6iiiiiiiiiiiiiiiiiixxxxhS xMMyMh xxhhhyMh xxxxx(ii).建立参数iM的方程组对 S(x)求导可得2211111( )() /22,6iiiiiiiiiiiiiixxxxSxMMyyhhhMMhxxx上

11、式中令ixx得 S(x) 在 xi 处的左导数()iSx，令1ixx得到右导数()iSx，因为 S(x) 在内节点 xi 处一阶导数连续，所以()()iiSxSx，进一步推导可得112,1,2,3,.,1iiiiiiMMMd in其中1iiiihhh，111iiiiihhh，精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 13 页1111116()6 ,1,2,.,1iiiiiiiiiiiiyyyydf xx xinhhhh上式为三弯矩方程组， 因为三弯矩方程组只有n-1 个方程，不能确定 n+1个未知量 Mi，所以需要再增加两个方程，

12、由边界条件确定。第一种边界条件 : 此时已知( )( )fafb和 . 不妨取0( )Mfa，( )nMfb，这时三弯矩方程组化为：1121101112111222,3,.,22iiiiiininnnnMMdMMMMdinMMdM以 上 方 程 组 系 数 矩阵 式 严 格 三 对 角 占 优 矩 阵 ， 可 用 追赶 法 求 解 。 求 出(1,2,.,1)iMin后，代入 S(x) 可得三次样条插值函数的数学表达式。第 二 种 边 界 条 件 ： 已 知( )( )fafb和。 记0( )( )nyfayfb,， 则 有00()()nnSxySxy,所以：1011101011 , ,366

13、3nnnnnnnnyyhhyyhhMMyMMyhh即0102MMd12nnnMMd其中10001116( )6( )nnnnnnyydyhhyydyhh所以得到第二种边界条件下的三弯矩方程组：0101112,2,1,2,3,.,1,2iiiiiinnnMMdMMMdinMMd该方程组系数矩阵是严格三对角占优矩阵，可用追赶法求解， 具体追赶法的求解过程见数值分析教材。第三种边界条件： 周期型边界条件 . 已知( )yf x是以0nTbaxx为周期的周期函数，则由周期性可知，01101111,nnnnnyyyy MMMMhh，精选学习资料 - - - - - - - - - 名师归纳总结 - -

14、- - - - -第 6 页，共 13 页这时，可以将点nx看成内点，则方程组对i=n 也成立，既有112nnnnnnMMMd，也即112nnnnnMMMd，其中11116()nnnnnnyyyydhhhh于是三弯矩方程组化为1121111112,2,2,3, 4,.,1,2.niiiiiinnnnnMMMdMMMdinMMMd该方程组可用 LU分解法求解，具体追赶法的求解过程见数值分析教材。程序框图如下：输出开始输入 , ( )( )(1),2,3.,h kx kx kkn( )( ) /( ( )(1),2,3,.1kh kh kh kknxy( ,2)nsize x1( )6 ( (1)

15、( ) / (1)( ( )(1)/ ( ) /( ( )(1),2,3,.1d ky ky kh ky ky kh kh kh kkn输入边界条件类型m1m分别输入和T( )Mn(1)M(2,2)Azeros nn( ,)2,( ,1)(1),(1, )(2)1,2,3.3A k kA k kkA kkkkn(2,2)2A nn(2,1)bzeros n( ,1)(1) ,2,3,.3b kd kkn(1,1)(2)(2)*(1)bdM(2,1)(1)(1)*( )b nd nnMn用追赶法求解AMb3321111211()()( )()666(),6iiiiiiiiiiiiiiiiiixx

16、xxhxxS xMMyMhhhhxxyMxxxhF2m( , )Azeros n n( , )2,( ,1)( ),(1, )(1)2,3.2A k kA k kkA kkkkn(1,1)2A( ,1)bzeros n( ,1)( ) ,1,2,3,.b kd kkn分别输入 f(a)和f(b) 一阶导数和yn0yT(1)6*(2)(1) / (2)0) /(2)dyyhyh( )6*( ( )(1) /( ) /( )d nyny ny nh nh n(1,2)1A(2,1)(2)A(1,1)2A nn( , )2A n n( ,1)1A n n(1, )(1)A nnn用追赶法求解AMb3

17、mFT(1,1)Azeros nn( ,)2,( ,1)(1),(1, )(2)1,2,3.2A k kA k kkA kkkkn(1,1)bzeros n( ,1)(1) ,1,2,3,.1b kd kkn( )6*(2)( ) /(2)( ( )(1) /( ) /( ( )(2)d nyy nhy ny nh nh nh(1,1)2A nn(1,1)(2)An(1,1)( )A nnF用LU分解法求解AMb结束精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 13 页程序使用说明： 因为该程序需要计算三种不同边界条件的三弯矩方程组，

18、所以首先定义追赶法和 LU 分解法求解线性方程组的函数，ZhuiGanFa(A,b) 和LU_fenjieqiuxianxingfangcheng(A,b)，然后在确定了边界条件类型之后，构造关于M 的三弯矩方程 AM=b 的A和b矩阵，然后分别代入函数 ZhuiGanFa(A,b) 和LU_fenjieqiuxianxingfangcheng(A,b)便可求得 M ，然后根据3321111211()()( )()666(),6iiiiiiiiiiiiiiiiiixxxxhxxS xMMyMhhhhxxyMxxxh便可得到，在1iixxx上的三次样条插值函数( )S x ，进而得到整个区间上的

19、三次样条差值函数( )S x 。算例计算结果：141页的计算实习当为第一类边界条件时：当-1=X=-0.8 时，S =0.3390+0.6443*X+0.4631*X2+0.1193*X3 当-0.8=X=-0.6 时，S =0.7658+2.245*X+2.464*X2+0.9529*X3 当-0.6=X=-0.4 时，S =0.7372+2.102*X+2.225*X2+0.8204*X3 当-0.4=X=-0.2 时，S =1.543+8.146*X+17.34*X2+13.41*X3 当-0.2=X=0时，S =-54.47*X3-0.2698e-14*X-23.39*X2+1.000

20、 当0=X=0.2时，S =1.000+0.2698e-14*X-23.39*X2+54.47*X3 当0.2=X=0.4时，S =1.543-8.146*X+17.34*X2-13.41*X3 当0.4=X=0.6时，S =0.7372-2.102*X+2.225*X2-0.8204*X3 当0.6=X=0.8时，S =0.7658-2.245*X+2.464*X2-0.9529*X3 当0.8=X=1时，S =0.3390-0.6443*X+0.4631*X2-0.1193*X3 当为第二类边界条件时：当-1=X=-0.8 时，S =0.3175+0.5663*X+0.3695*X2+0.

21、8225e-1*X3 当-0.8=X=-0.6 时，S =0.7683+2.257*X+2.483*X2+0.9628*X3 当-0.6=X=-0.4 时，S =0.7370+2.100*X+2.222*X2+0.8177*X3 当-0.4=X=-0.2 时，S =1.543+8.146*X+17.34*X2+13.41*X3 当-0.2=X=0时，S =-54.47*X3-0.2320e-14*X-23.39*X2+1.000 -1-0.8-0.6-0.4-0.200.20.40.60.8100.10.20.30.40.50.60.70.80.91xy三 次 样 条 差 值 函 数 S10(

22、x) 和 原 函 数 f(x) 的 对 比红 线 ： 原 函 数蓝 线 ： S10(x)精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 13 页当0=X=0.2时，S =1.000+0.2043e-14*X-23.39*X2+54.47*X3 当0.2=X=0.4时，S =1.543-8.146*X+17.34*X2-13.41*X3 当0.4=X=0.6时，S =0.7370-2.100*X+2.222*X2-0.8177*X3 当0.6=X=0.8时,S =0.7683-2.257*X+2.483*X2-0.9628*X3 当0.8

23、=X=1时,S =0.3175-0.5663*X+0.3695*X2-0.8225e-1*X3 3. 龙贝格积分法：对于复化梯形求积公式，取积分近似值2221()()4 1nnnnI fTTTT对复化辛普森求积公式4(4)()( ),2880nbaI fSh fab4(4)211()( )(),2880 2nba hI fSfab因为(4)(4)1( )()ff，所以上述两式相除得2()16()nnIfSIfS所以，22221( )()41nnnnI fSSSS同理，22231( )()41nnnnI fCCCC对2nT，2nS和2nC分析可得222222231()4 11()411()41n

24、nnnnnnnnnnnSTTTCSSSRCCC-1-0.8-0.6-0.4-0.200.20.40.60.8100.10.20.30.40.50.60.70.80.91xy红 线 ： 原 函 数蓝 线 ： S10(x)三 次 样 条 差 值 函 数 S10(x) 和 原 函 数 f(x) 的 对 比精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 13 页龙贝格积分算法如下：1111111121122122222222232222( )( ),21(21),0,1,2.,2221(),4 11(),0,1,2.,411(),41kkkkk

25、kkkkkkkkkkkkibaTf af bbabaTTf aikSTTTCSSSkRCCC如果122,kkRR则取12kI fR，否则，继续计算直到满足条件。程序框图：开始输入 , , , b1eps3k1kkfa23221kkRRepsT输出近似解22kIfRF结束1( )( )2baTf af b12112211(21),0,1,2222kkkkkibabaTTf aik1122221(),0,1,241kkkkSTTTk11222221(),0,141kkkkCSSSk11322221() ,041kkkkRCCCk12112211(21)222kkkkkibabaTTf ai1122

26、221()41kkkkSTTT11222221()41kkkkCSSS2112322221()41kkkkRCCC精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 13 页程序使用说明 : 运行本程序的时候，直接按照提示输入所求积分的原函数( )f x（比如11x） ，然后按提示依次输入积分下限a，积分上限 b 和积分精度1eps，然后程序便可计算出原函数( )f x在 , a b之间的积分数值。算例计算结果：6.2 1010.69266061dxx120ln(1)0.27291381xdxx10ln(1)0.8222677xdxx/2

27、0sin( )1.3714136xdxx4. 四阶龙格 -库塔法求解常微分方程的初值问题算法原理：用标准 4级4阶R-K法求解一阶常微分方程的算法公式为：1123412132431(22)6(,)11(,)2211(,)22(,)iiiiiiiiiiyyKKKKKhf x yKhf xh yKKhf xh yKKhf xh yK求解m 阶常微分方程()(1)(1)(1)000( , ,.,);,( ),( ) ,.,( )mmmmyf x y y yyaxby ayy ayyay将其转化为一阶微分方程组来求解。 为此， 引进新的变量1yy , 2yy, , (1)mmyy, 即可将 m 阶微分

28、方程，转化为如下的一阶微分方程组：1223112,.,(1)10200.( ,)( ),( ) ,.,( )mmmmmmyyyyyyyf x y yyy ayy ayyay程序框图：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 13 页开始输入 , , , bhma结束1m()/mbahxa11(1)(1,1)yy输入(1,1)yT1(1,1)K(1,1)=h*(eval(ym)x=x+h/2y(1,1)=y1+K(1,1)/2y1=y(1,1)K(1,2)=h*(eval(ym) y(1,1)=y1+K(1,2)/2-K(1,1)

29、/2y1=y(1,1)K(1,3)=h*(eval(ym)x=x+h/2y(1,1)=y1+K(1,3)-K(1,2)/2y1=y(1,1)K(1,4)=h*(eval(ym)y11yy(k1)=y11(k1-1)+(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4)/6y(1,1)=y11(k1)x=a+(k1-1)*h12k11kmm11 1kkF输出T11y输入y(k2,1)21:km() /y22(1,k2)=y(k2,1)k2=1:mK(k,1)=h*y(k+1,1)1:(1)K(m,1)=h*(eval(ym)x=x+h/2y(k3,1)=y(k3,1)+K(k3,1)

30、/231:K(k,2)=h*y(k+1,1)k=1:(m-1)K(m,2)=h*(eval(ym)y(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2k3mbahxakmkm=1:mK(k,3)=h*y(k+1,1)k=1:(m-1)K(m,3)=h*(eval(ym)x=x+h/2y(k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2k3=1:mK(k,4)=h*y(k+1,1)k=1:(m-1)K(m,4)=h*(eval(ym)y22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4)/6k5=1

31、:my(k6,1)=y22(k4,k6)k6=1:mx=a+(k4-1)*h441kk42kF41kmmF输出y22(:,1)F程序使用说明：运行该程序，按照提示依次输入常微分方程的阶数m ，x的下限a和上限 b，步长 h以及y的最高阶导数表达式()(1)(2).( )mmmyaybyf x， 特别需要注意的是 ： 由于程序设定的问题，一定要将公式中的()ky用(1,1)y k代替，例如所求解的初值问题为2 22sin01xyyyexx，那么本程序所需要输入的公式为：2*(2,1)2*(1,1) exp(2*)*sin()yyxx。且( )f x代表含x的表达式。接着依次输入已知条件( )y

32、a，( )y a ，( )ya ，即可，最终可以得出该问题的精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 13 页解。算例计算结果：9.2 题结果(0)-1.000000000000000(0.05)-0.847436458333333(0.1)-0.689482936121419(0.15)-0.525724908067489(0.2)-0.355719511320969(0.25)-0.178993729355764(0.3)0.004957535888930(0.35)0.196673510288970yyyyyyyy(0.4

33、)0.396729719312157(0.45)0.605740292781225(0.5)0.824360410550626(0.55)1.053288897488131(0.6)1.293270976642268(0.65)1.545101189994835(0.7)1.809626496745704(0.75)2.087749559656611(0.yyyyyyyyy8)2.380432230593372(0.85)2.688699247053874(0.9)3.013642152154253(0.95)3.356423451269989(1)3.718281019294475yyyy00.10.20.30.40.50.60.70.80.91-1-0.500.511.522.533.54xy红 色 +： R-K近 似 解蓝 线 ： 解 析 解近 似 解 和 解 析 解 y=x*ex+2*x-1 的 比 较精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 13 页