西交计算方法A上机大作业

计算方法A 上机大作业1. 共轭梯度法求解线性方程组算法原理：由定理3.4.1可知系数矩阵A 是对称正定矩阵的线性方程组Ax=b 的解与求解二次函数1()2TT f x x Ax b x =- 极小点具有等价性，所以可以利用共轭梯度法求解1()2TT f x x Ax b x =-的极小点来达到求解Ax=b 的目的。

共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式：(1)()()k k k k x x d α+=+产生的迭代序列(1)(2)(3)x x x ，，，... 在无舍入误差假定下，最多经过n 次迭代，就可求得()f x 的最小值，也就是方程Ax=b 的解。

首先导出最佳步长k α的计算式。

假设迭代点()k x 和搜索方向()k d 已经给定，便可以通过()()()()k k f x d φαα=+的极小化()()min ()()k k f x d φαα=+来求得，根据多元复合函数的求导法则得:()()()'()()k k T k f x d d φαα=∇+令'()0φα=，得到:()()()()k T k k k T k r d d Adα= ，其中()()k k r b Ax =-然后确定搜索方向()k d 。

给定初始向量(0)x 后，由于负梯度方向是函数下降最快的方向，故第一次迭代取搜索方向(0)(0)(0)(0)()dr f x b Ax ==-∇=- 。

令(1)(0)00x x d α=+其中(0)(0)0(0)(0)T T r d d Adα=。

第二次迭代时，从(1)x 出发的搜索方向不再取(1)r ，而是选取(1)(1)(0)0dr d β=+，使得(1)d 与(0)d 是关于矩阵A 的共轭向量，由此可求得参数0β:(1)(0)0(0)(0)T T r Ad d Adβ=-然后从(1)x 出发，沿(1)d 进行搜索得到(2)(1)(1)1x x d α=+设已经求出(1)()()k k k k x x d α+=+，计算(1)(1)k k r b Ax ++=-。

令(1)(1)()k k k k dr d β++=+，选取k β，使得(1)k d +和()k d 是关于A 的共轭向量，可得:(1)()()()k T k k k T k r Ad d Adβ+=-具体编程计算过程如下：（1） 给定初始近似向量(0)x 以及精度0ε> ； （2） 计算(0)(0)r b Ax =-，取(0)(0)d r =； （3） For k=0 to n-1 do（i ）()()0()()k T k k T k r d d Adα=；（ii ）(1)()()k k k k xx d α+=+；（iii ）(1)(1)k k r b Ax ++=-；（iv ）若(1)k r ε+≤或k=n-1，则输出近似解(1)k x + ，停止；否则，转（v ）；（v ）2(1)22()2k k k r rβ+= ；（vi ）(1)(1)()k k k k d r d β++=+；End do程序框图：程序使用说明：本共轭梯度法求解线性方程的程序直接打开matlab运行，在求解线性方程组Ax=b（A是对称正定矩阵）的时候，直接运行程序Gongetidufa,输入A,b的值，虽然该函数是通用的，但是对于矩阵A和向量b的输入，需要使用者根据A和b的特点自行输入。

算例3.4.2计算结果：对99页例题3.4.2，运行程序Gongetidufa将矩阵A,b读入系统程序如下：clear allclcA=input('请输入A的值'); %输入n阶正定矩阵A的值b=input('请输入b的值:'); %输入b的值n=size(A,1); %求出矩阵A的行数x=zeros(n,1); %给定x的初始值e=10^(-12); %给出精度r=b-A*x;d=r;for(i=1:n)a=norm(r,2)^2/(d'*A*d); %求最佳步长 x=x+a*d; j=r;r=b-A*x;if(norm(r)<=e||i==n) break; elseB=norm(r,2)^2/norm(j,2)^2; d=r+B*d; end endx %输出最终的x 的结果 计算结果：x=[1;1;1]2.三次样条差值算法原理（三次样条插值函数的导出）： (i).导出在子区间[]1,i i x x - 上的S(x)的表达式由于S(x)的二阶导数连续，设S(x)再节点i x 处的二阶导数值S ’’(xi)=Mi ，其中Mi 为未知的待定参数。

由S(x)是分段的三次多项式知，S ’’(x)是分段线性函数，S ’’(x)在子区间[]1,i i x x -上可表示为1111111''(),i i i ii i i i i i i i i ii ix x x x S x M M x x x x x x x x M M x x x h h ---------=+----=+≤≤其中hi=xi-x(i-1),对上式两次积分得到()()331111()66()(),i i i i iii i i i i ix x x x S x M M h h b x x c x x x x x ------=++-+-≤≤由插值条件11(),()i i i i S x y S x y --== 得到2211()/,()/66i i i i i i i i i i h h b y M h c y M h --=-=-将i i b c 和 代入()S x 可得()()3321111211()()/()666()/(),6ii i i i i i i i i iii i i i i i x x x x h S x M M y M h x x h h hy M h x x x x x --------=++--+--≤≤(ii).建立参数i M 的方程组 对S(x)求导可得()()2211111'()()/22,6ii i i i i iiii i i i ix x x x S x M M y y h h h M M h x x x -------=-++---≤≤上式中令i x x = 得S(x)在xi 处的左导数'()i S x - ，令1i x x -=得到右导数'()i S x +，因为S(x)在内节点xi 处一阶导数连续，所以''()()i i S x S x -+=，进一步推导可得112,1,2,3,...,1i i i i i i M M M d i n μλ-+++==-其中1ii i i h h h μ+=+，111i i i i i h h h λμ++==-+，1111116()6[,,],1,2,...,1i i i i i i i i i i i iy y y y d f x x x i n h h h h +--+++--=-==-+上式为三弯矩方程组，因为三弯矩方程组只有n-1个方程，不能确定n+1个未知量Mi ，所以需要再增加两个方程，由边界条件确定。

第一种边界条件:此时已知''()''()f a f b 和 .不妨取0''()M f a =，''()n M f b =，这时三弯矩方程组化为：1121101112111222,3,...,22i i i i i in i n n n n M M d M M M M d i n M M d Mλμμλμλ-+-----+=-⎧⎪++==-⎨⎪+=-⎩ 以上方程组系数矩阵式严格三对角占优矩阵，可用追赶法求解。

求出(1,2,...,1)i M i n =- 后，代入S(x)可得三次样条插值函数的数学表达式。

第二种边界条件：已知'()'()f a f b 和。

记0''()''()n y f a y f b ==,，则有00'()''()'n n S x y S x y +-==,所以：1011101011',',3663n n n n n n n ny y h h y y h hM M y M M y h h ------+=++= 即0102M M d +=12n n n M M d -+=其中10001116(')6(')n n n n n ny y d y h h y y d y h h --=--=-所以得到第二种边界条件下的三弯矩方程组：0101112,2,1,2,3,...,1,2i i i i i i n n n M M d M M M d i n M M dμλ-+-+=⎧⎪++==-⎨⎪+=⎩ 该方程组系数矩阵是严格三对角占优矩阵，可用追赶法求解，具体追赶法的求解过程见《数值分析》教材。

第三种边界条件：周期型边界条件.已知()y f x = 是以0n T b a x x =-=- 为周期的周期函数，则由周期性可知，01101111,,,,n n n n n y y y y M M M M h h +++===== ，这时，可以将点n x 看成内点，则方程组对i=n 也成立，既有112n n n n n n M M M d μλ-+++= ，也即112n n n n n M M M d λμ-++=，其中11116()n n n n n ny y y y d h h h h ---=-+ 于是三弯矩方程组化为1121111112,2,2,3,4,....,1,2.n i i i i i i n n n n n M M M d M M M d i n M M M d λμμλλμ-+-++=⎧⎪++==-⎨⎪++=⎩ 该方程组可用matlab 直接解出。

程序框图如下：程序使用说明：本程序是求解137页例题4.6.1的运行结果，通过程序便可求得M ，然后根据3321111211()()()()666(),6i i i ii i i i i i iii i i i i ix x x x h x x S x M M y M h h h h x x y M x x x h ---------=++--+-≤≤便可得到，在1i i x x x -≤≤ 上的三次样条插值函数()S x ，进而得到整个区间上的三次样条差值函数()S x 。

算例计算结果：137页例题4.6.1的计算实习1、打开matlab 运行Sanciyangtiao 程序2、自行输入x 和y 的节点值3、得出计算结果3.龙贝格积分法对于复化梯形求积公式，取积分近似值2221()()41n n n n I f T T T T ≈+-=- 对复化辛普森求积公式4(4)()(),2880n b a I f S h f a b ηη--≈-≤≤4(4)211()()(),28802n b a h I f S f a b ηη--≈-≤≤ 因为(4)(4)1()()ff ηη=，所以上述两式相除得2()16()nnI f S I f S -≈-所以，22221()()41n n n n I f S S S S ≈+-=-同理，22231()()41n n n n I f C C C C ≈+-=-对2n T ，2n S 和2n C 分析可得222222231()411()411()41n n n n n n n n n n n n S T T T C S S S R C C C ⎧=+-⎪-⎪⎪=+-⎨-⎪⎪=+-⎪-⎩龙贝格积分算法如下：1111111121122122222222232222[()()],21((21)),0,1,2...,2221(),411(),0,1,2...,411(),41kk k k k k k k k k k k k k k k k i b a T f a f b b a b aT T f a i k S T T T C S S S k R C C C +++++++++=-=+--=++-=⎧=+-⎪-⎪⎪=+-=⎨-⎪⎪=+-⎪-⎩∑如果122,k k R R ε+-< 则取12[]k I f R +≈，否则，继续计算直到满足条件。