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

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: 精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 1 页，共 13 页 - - - - - - - - -精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 1 页，共 13 页 - - - - - - - - -(1)(0)0(0)(0)ttraddad然后从(1)x出发，沿(1)d进行搜索得到(2)(1)(1)1xxd设已经求出(1)( )( )kkkkxxd，计算(1)(1)kkrbax

4、。令(1)(1)( )kkkkdrd，选取k，使得(1)kd和( )kd是关于 a 的共轭向量，可得: (1)()( )( )ktkkk tkraddad具体编程计算过程如下：（1） 给定初始近似向量(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

5、)()kkkkdrd；end do 程序框图：精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 2 页，共 13 页 - - - - - - - - -精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 2 页，共 13 页 - - - - - - - - -开始输入 , , , (0)x(0)(0)rbax(0)(0)dr( )( )( )( )k tkkk tkrddad(1)( )()kkkkxxd(1)(1)kkrbax0k1kkab( ,1)nsize a或(1)kr1knt输出近似解(1)kxf结束2

6、(1)22()2kkkrr(1)(1)( )kkkkdrd程 序 使 用 说 明 ： 本 共 轭 梯 度 法 求 解 线 性 方 程 的 程 序 是 一 个 函 数gongetidu2(a,b,x0,epsa)，在求解线性方程组ax=b（a 是对称正定矩阵）的时候 ， 首 先 给 定 初 始 向 量x0和 误 差epsa ， 然 后 直 接 调 用 该 函 数gongetidu2(a,b,x0,epsa)即可，其中 a,b,x0 和 epsa 都是可变的，虽然该函数是通用的，但是对于矩阵 a和向量 b 的输入，需要使用者根据 a和 b 的特点自行输入。算例 3.2 计算结果：对题 113 页

7、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 的矩阵 a a(k,k)=-2; a(k,k+1)=1; a(k+1,k)=1; end 精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 3 页，共 13 页 - - - - - - - -

8、-精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 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=gongetidu2(a,b,x0,epsa) n=size(a,1

9、); 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;.1;1;1;1;1;1;1;1;1;1 2. 三次

10、样条差值算法原理（三次样条插值函数的导出） ：(i).导出在子区间1,iixx上的 s(x) 的表达式由于 s(x) 的二阶导数连续，设 s(x) 再节点ix处的二阶导数值 s(xi)=mi ，其中 mi 为未知的待定参数。由s(x) 是分段的三次多项式知，s(x) 是分段线性函数， s(x) 在子区间1,iixx上可表示为精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 4 页，共 13 页 - - - - - - - - -精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 4 页，共 13 页 - - -

11、- - - - - -1111111( ),iiiiiiiiiiiiiiiixxxxsxmmxxxxxxxxmmxxxhh其中 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的方

12、程组对 s(x)求导可得2211111( )() /22,6iiiiiiiiiiiiiixxxxsxmmyyhhhmmhxxx上式中令ixx得 s(x) 在 xi 处的左导数()isx，令1ixx得到右导数()isx，因为 s(x) 在内节点 xi 处一阶导数连续，所以()()iisxsx，进一步推导可得112,1,2,3,.,1iiiiiimmmd in其中1iiiihhh，111iiiiihhh，精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 5 页，共 13 页 - - - - - - - - -精品学习资料 可选择p d f - - - -

13、- - - - - - - - - - 第 5 页，共 13 页 - - - - - - - - -1111116()6 ,1,2,.,1iiiiiiiiiiiiyyyydf xx xinhhhh上式为三弯矩方程组， 因为三弯矩方程组只有n-1 个方程，不能确定 n+1个未知量 mi，所以需要再增加两个方程，由边界条件确定。第一种边界条件 : 此时已知( )( )fafb和 . 不妨取0( )mfa，( )nmfb，这时三弯矩方程组化为：1121101112111222,3,.,22iiiiiininnnnmmdmmmmdinmmdm以 上 方 程 组 系 数 矩阵 式 严 格 三 对 角 占

14、 优 矩 阵 ， 可 用 追赶 法 求 解 。 求 出(1,2,.,1)imin后，代入 s(x) 可得三次样条插值函数的数学表达式。第 二 种 边 界 条 件 ： 已 知( )( )fafb和。 记0( )( )nyfayfb,， 则 有00()()nnsxysxy,所以：1011101011 , ,3663nnnnnnnnyyhhyyhhmmymmyhh即0102mmd12nnnmmd其中10001116( )6( )nnnnnnyydyhhyydyhh所以得到第二种边界条件下的三弯矩方程组：0101112,2,1,2,3,.,1,2iiiiiinnnmmdmmmdinmmd该方程组系数矩

15、阵是严格三对角占优矩阵，可用追赶法求解， 具体追赶法的求解过程见数值分析教材。第三种边界条件： 周期型边界条件 . 已知( )yf x是以0ntbaxx为周期的周期函数，则由周期性可知，01101111,nnnnnyyyy mmmmhh，精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 6 页，共 13 页 - - - - - - - - -精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 6 页，共 13 页 - - - - - - - - -这时，可以将点nx看成内点，则方程组对i=n 也成立，既有112n

16、nnnnnmmmd，也即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)( ) / (1)( ( )(1)/ ( ) /( ( )(1),2,3,.1d ky ky kh ky k

17、y 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(),6iiiiiiiiiiiiiiiiiixxxxhxxs xmmymhhhhxxymxxxhf2m( , )azeros n n( , )2,( ,1)

18、( ),(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用追赶法求解amb3mft(1,1)azeros nn( ,)2,( ,1)(1),(1, )(2)1,2,3.2a k ka

19、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结束精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 7 页，共 13 页 - - - - - - - - -精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 7 页，共 13 页 - - - - - - -

20、 - -程序使用说明： 因为该程序需要计算三种不同边界条件的三弯矩方程组，所以首先定义追赶法和 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便可得

21、到，在1iixxx上的三次样条插值函数( )s x ，进而得到整个区间上的三次样条差值函数( )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 =

22、-54.47*x3-0.2698e-14*x-23.39*x2+1.000 当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=-

23、0.8 时，s =0.3175+0.5663*x+0.3695*x2+0.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

24、.50.60.70.80.91xy三 次 样 条 差 值 函 数 s10(x) 和 原 函 数 f(x) 的 对 比红 线 ： 原 函 数蓝 线 ： s10(x)精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 8 页，共 13 页 - - - - - - - - -精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 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-

25、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=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)(

26、4)1( )()ff，所以上述两式相除得2()16()nnifsifs所以，22221( )()41nnnni fssss同理，22231( )()41nnnni fcccc对2nt，2ns和2nc分析可得222222231()4 11()411()41nnnnnnnnnnnnstttcsssrccc-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) 的 对 比精品学习资料 可选择p d f - - -

27、- - - - - - - - - - - 第 9 页，共 13 页 - - - - - - - - -精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 9 页，共 13 页 - - - - - - - - -龙贝格积分算法如下：1111111121122122222222232222( )( ),21(21),0,1,2.,2221(),4 11(),0,1,2.,411(),41kkkkkkkkkkkkkkkkkibatf af bbabattf aikstttcssskrccc如果122,kkrr则取12ki fr，否则，继续计算直到满足条件。程

28、序框图：开始输入 , , , b1eps3k1kkfa23221kkrrepst输出近似解22kifrf结束1( )( )2batf af b12112211(21),0,1,2222kkkkkibabattf aik1122221(),0,1,241kkkkstttk11222221(),0,141kkkkcsssk11322221() ,041kkkkrccck12112211(21)222kkkkkibabattf ai1122221()41kkkksttt11222221()41kkkkcsss2112322221()41kkkkrccc精品学习资料 可选择p d f - - - -

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

30、22677xdxx/20sin( )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

31、, 即可将 m 阶微分方程，转化为如下的一阶微分方程组：1223112,.,(1)10200.( ,)( ),( ) ,.,( )mmmmmmyyyyyyyf x y yyy ayy ayyay程序框图：精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 11 页，共 13 页 - - - - - - - - -精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 11 页，共 13 页 - - - - - - - - -开始输入 , , , bhma结束1m()/mbahxa11(1)(1,1)yy输入(1,1)y

32、t1(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)/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

33、)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)/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:my(k6,1)=y22(k4,k6)k6=1:mx=a+(k4-1)*h441kk42kf41kmmf输出y22(:,1)f程序使用说明：运行该程序，按照提示依次输入常微分方程的阶数m ，x的下限a和上限 b，步长 h以及y的最高阶导数表达式()(1)