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

计算方法A上机大作业1.共轭梯度法求解线性方程组算法原理：由定理3.4.1可知系数矩阵A是对称正定矩阵的线性方程组Ax=b的解与求解二次函数1()2TTfxxAxbx极小点具有等价性，所以可以利用共轭梯度法求解1()2TTfxxAxbx的极小点来达到求解Ax=b的目的。共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式：(1)()()kkkkxxd产生的迭代序列(1)(2)(3)xxx，，，...在无舍入误差假定下，最多经过n次迭代，就可求得()fx的最小值，也就是方程Ax=b的解。首先导出最佳步长k的计算式。假设迭代点()kx和搜索方向()kd已经给定，便可以通过()()()()kkfxd的极小化()()min()()kkfxd来求得，根据多元复合函数的求导法则得:()()()'()()kkTkfxdd令'()0，得到:()()()()kTkkkTkrddAd，其中()()kkrbAx然后确定搜索方向()kd。给定初始向量(0)x后，由于负梯度方向是函数下降最快的方向，故第一次迭代取搜索方向(0)(0)(0)(0)()drfxbAx。令(1)(0)00xxd其中(0)(0)0(0)(0)TTrddAd。第二次迭代时，从(1)x出发的搜索方向不再取(1)r，而是选取(1)(1)(0)0drd，使得(1)d与(0)d是关于矩阵A的共轭向量，由此可求得参数0:(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)()()()kTkkkTkrAddAd具体编程计算过程如下：（1）给定初始近似向量(0)x以及精度0；（2）计算(0)(0)rbAx，取(0)(0)dr；（3）Fork=0ton-1do（i）()()0()()kTkkTkrddAd；（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；Enddo程序框图：开始输入,,,(0)x(0)(0)rbAx(0)(0)dr()()()()kTkkkTkrddAd(1)()()kkkkxxd(1)(1)kkrbAx0k1kkAb(,1)nsizeA或(1)kr1knT输出近似解(1)kxF结束2(1)22()2kkkrr(1)(1)()kkkkdrd程序使用说明：本共轭梯度法求解线性方程的程序直接打开matlab运行，在求解线性方程组Ax=b（A是对称正定矩阵）的时候，直接运行程序Gongetidufa,输入A,b的值，虽然该函数是通用的，但是对于矩阵A和向量b的输入，需要使用者根据A和b的特点自行输入。算例3.4.2计算结果：对99页例题3.4.2，运行程序Gongetidufa将矩阵A,b读入系统程序如下：clearallclcA=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;endendx%输出最终的x的结果计算结果：x=[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上可表示为1111111''(),iiiiiiiiiiiiiiiixxxxSxMMxxxxxxxxMMxxxhh其中hi=xi-x(i-1),对上式两次积分得到331111()66()(),iiiiiiiiiiiixxxxSxMMhhbxxcxxxxx由插值条件11(),()iiiiSxySxy得到2211()/,()/66iiiiiiiiiihhbyMhcyMh将iibc和代入()Sx可得3321111211()()/()666()/(),6iiiiiiiiiiiiiiiiiixxxxhSxMMyMhxxhhhyMhxxxxx(ii).建立参数iM的方程组对S(x)求导可得2211111'()()/22,6iiiiiiiiiiiiiixxxxSxMMyyhhhMMhxxx上式中令ixx得S(x)在xi处的左导数'()iSx，令1ixx得到右导数'()iSx，因为S(x)在内节点xi处一阶导数连续，所以''()()iiSxSx，进一步推导可得112,1,2,3,...,1iiiiiiMMMdin其中1iiiihhh，111iiiiihhh，1111116()6[,,],1,2,...,1iiiiiiiiiiiiyyyydfxxxinhhhh上式为三弯矩方程组，因为三弯矩方程组只有n-1个方程，不能确定n+1个未知量Mi，所以需要再增加两个方程，由边界条件确定。第一种边界条件:此时已知''()''()fafb和.不妨取0''()Mfa，''()nMfb，这时三弯矩方程组化为：1121101112111222,3,...,22iiiiiininnnnMMdMMMMdinMMdM以上方程组系数矩阵式严格三对角占优矩阵，可用追赶法求解。求出(1,2,...,1)iMin后，代入S(x)可得三次样条插值函数的数学表达式。第二种边界条件：已知'()'()fafb和。记0''()''()nyfayfb,，则有00'()''()'nnSxySxy,所以：1011101011',',3663nnnnnnnnyyhhyyhhMMyMMyhh即0102MMd12nnnMMd其中10001116(')6(')nnnnnnyydyhhyydyhh所以得到第二种边界条件下的三弯矩方程组：0101112,2,1,2,3,...,1,2iiiiiinnnMMdMMMdinMMd该方程组系数矩阵是严格三对角占优矩阵，可用追赶法求解，具体追赶法的求解过程见《数值分析》教材。第三种边界条件：周期型边界条件.已知()yfx是以0nTbaxx为周期的周期函数，则由周期性可知，01101111,,,,nnnnnyyyyMMMMhh，这时，可以将点nx看成内点，则方程组对i=n也成立，既有112nnnnnnMMMd，也即112nnnnnMMMd，其中11116()nnnnnnyyyydhhhh于是三弯矩方程组化为1121111112,2,2,3,4,....,1,2.niiiiiinnnnnMMMdMMMdinMMMd该方程组可用matlab直接解出。程序框图如下：输出开始输入,()()(1),2,3...,hkxkxkkn()()/(()(1)),2,3,...1khkhkhkknxy(,2)nsizex1()6(((1)())/(1)(()(1))/())/(()(1)),2,3,...1dkykykhkykykhkhkhkkn输入边界条件类型m1m分别输入和T()Mn(1)M(2,2)Azerosnn(,)2,(,1)(1),(1,)(2)1,2,3...3AkkAkkkAkkkkn(2,2)2Ann(2,1)bzerosn(,1)(1),2,3,...3bkdkkn(1,1)(2)(2)*(1)bdM(2,1)(1)(1)*()bndnnMn用追赶法求解AMb3321111211()()()()666(),6iiiiiiiiiiiiiiiiiixxxxhxxSxMMyMhhhhxxyMxxxhF2m(,)Azerosnn(,)2,(,1)(),(1,)(1)2,3...2AkkAkkkAkkkkn(1,1)2A(,1)bzerosn(,1)(),1,2,3,...bkdkkn分别输入f(a)和f(b)一阶导数和yn0yT(1)6*(((2)(1))/(2)0)/(2)dyyhyh()6*((()(1))/())/()dnynynynhnhn(1,2)1A(2,1)(2)A(1,1)2Ann(,)2Ann(,1)1Ann(1,)(1)Annn用追赶法求解AMb3mFT(1,1)Azerosnn(,)2,(,1)(1),(1,)(2)1,2,3...2AkkAkkkAkkkkn(1,1)bzerosn(,1)(1),1,2,3,...1bkdkkn()6*(((2)())/(2)(()(1))/())/(()(2))dnyynhynynhnhnh(1,1)2Ann(1,1)(2)An(1,1)()AnnF用LU分解法求解AMb结束程序使用说明：本程序是求解137页例题4.6.1的运行结果，通过程序便可求得M，然后根据3321111211()()()()666(),6iiiiiiiiiiiiiiiiiixxxxhxxSxMMyMhhhhxxyMxxxh便可得到，在1iixxx上的三次样条插值函数()Sx，进而得到整个区间上的三次样条差值函数()Sx。算例计算结果：137页例题4.6.1的计算实习1、打开matlab运行Sanciyangtiao程序2、自行输入x和y的节点值3、得出计算结果3.龙贝格积分法对于复化梯形求积公式，取积分近似值2221()()41nnnnIfTTTT对复化辛普森求积公式4(4)()(),2880nbaIfShfab4(4)211()()(),28802nbahIfSfab因为(4)(4)1()()ff，所以上述两式相除得2()16()nnIfSIfS所以，22221()()41nnnnIfSSSS同理，22231()()41nnnnIfCCCC对2nT，2nS和2nC分析可得222222231()411()411()41nnnnnnnnnnnnSTTTCSSSRCCC龙贝格积分算法如下：1111111121122122222222232222[