《2022年西安交通大学2021年计算方法A上机大作业》由会员分享,可在线阅读,更多相关《2022年西安交通大学2021年计算方法A上机大作业(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)00xxd其中(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 , 然 后 直 接 调 用 该 函 数Go
6、ngetidu2(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=Gongetid
8、u2(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其中 h
10、i=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 , ,3663
13、nnnnnnnnyyhhyyhhMMyMMyhh即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(),6iiiiiiiiiiiiiiiiiixxx
16、xhxxS 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用追赶法求解AMb3m
17、FT(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.8
21、225e-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(x
22、) 和 原 函 数 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()41nn
24、nnnnnnnnnnSTTTCSSSRCCC-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(),41kkkkkk
25、kkkkkkkkkkkibaTf 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 ai11222
26、21()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/20
27、sin( )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 a
32、,( )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 页