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

计算方法A上机大作业张晓璐 硕4011班 学号：31140090971. 共轭梯度法求解线性方程组算法原理：由定理3.4.1可知系数矩阵A是对称正定矩阵的线性方程组Ax=b的解与求解二次函数 极小点具有等价性，所以可以利用共轭梯度法求解的极小点来达到求解Ax=b的目的。共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式：产生的迭代序列 在无舍入误差假定下，最多经过n次迭代，就可求得 的最小值，也就是方程Ax=b的解。首先导出最佳步长的计算式。假设迭代点和搜索方向已经给定，便可以通过 的极小化来求得，根据多元复合函数的求导法则得:令，得到: ，其中然后确定搜索方向。给定初始向量后，由于负梯度方向是函数下降最快的方向，故第一次迭代取搜索方向 。令其中。第二次迭代时，从 出发的搜索方向不再取，而是选取，使得与是关于矩阵A的共轭向量，由此可求得参数:然后从出发，沿进行搜索得到设已经求出，计算。令，选取，使得和是关于A的共轭向量，可得:具体编程计算过程如下：（1） 给定初始近似向量 以及精度 ；（2） 计算，取；（3） For k=0 to n-1 do（i）；（ii）；（iii）；（iv）若，则输出近似解 ，停止；否则，转（v）； （v） ；（vi）；End do程序框图：程序使用说明：本共轭梯度法求解线性方程的程序是一个函数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页3.2，首先编程将矩阵A,b,x0,epsa读入系统，然后再调用函数Gongetidu2(A,b,x0,epsa)；程序如下：clear A b x0 %程序运行前首先清除A，b和X0的数值，以免影响计算clcn=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;endA(n,n)=-2;A;b(1,1)=-1; %读取题目3.2的矩阵bb(n,1)=-1;b;x0(1,1)=input(假定初始向量各元素相同，均为：); %给题目3.2的迭代过程赋初值for kk=2:n x0(kk,1)=x0(1,1);endx=Gongetidu2(A,b,x0,epsa) %调用共轭梯度法求线性方程组的函数 function x=Gongetidu2(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;.1;1;1;1;1;1;1;1;1;12. 三次样条差值算法原理（三次样条插值函数的导出）：(i).导出在子区间 上的S(x)的表达式由于S(x)的二阶导数连续，设S(x)再节点处的二阶导数值S(xi)=Mi，其中Mi为未知的待定参数。由S(x)是分段的三次多项式知，S(x)是分段线性函数，S(x)在子区间上可表示为其中hi=xi-x(i-1),对上式两次积分得到由插值条件 得到将 代入 可得(ii).建立参数 的方程组对 S(x)求导可得上式中令 得S(x)在xi处的左导数 ，令得到右导数，因为S(x)在内节点xi处一阶导数连续，所以，进一步推导可得 其中，上式为三弯矩方程组，因为三弯矩方程组只有n-1个方程，不能确定n+1个未知量Mi，所以需要再增加两个方程，由边界条件确定。第一种边界条件:此时已知 .不妨取，这时三弯矩方程组化为：以上方程组系数矩阵式严格三对角占优矩阵，可用追赶法求解。求出 后，代入S(x)可得三次样条插值函数的数学表达式。第二种边界条件：已知。记，则有所以：即其中所以得到第二种边界条件下的三弯矩方程组：该方程组系数矩阵是严格三对角占优矩阵，可用追赶法求解，具体追赶法的求解过程见数值分析教材。第三种边界条件：周期型边界条件.已知 是以 为周期的周期函数，则由周期性可知， ，这时，可以将点 看成内点，则方程组对i=n也成立，既有 ，也即，其中于是三弯矩方程组化为该方程组可用LU分解法求解，具体追赶法的求解过程见数值分析教材。 程序框图如下：程序使用说明：因为该程序需要计算三种不同边界条件的三弯矩方程组，所以首先定义追赶法和LU分解法求解线性方程组的函数，ZhuiGanFa(A,b)和LU_fenjieqiuxianxingfangcheng(A,b)，然后在确定了边界条件类型之后，构造关于M的三弯矩方程AM=b的A和b矩阵，然后分别代入函数ZhuiGanFa(A,b)和LU_fenjieqiuxianxingfangcheng(A,b)便可求得M，然后根据便可得到，在 上的三次样条插值函数，进而得到整个区间上的三次样条差值函数。算例计算结果：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 当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.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当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=X=1时,S =0.3175-0.5663*X+0.3695*X2-0.8225e-1*X33. 龙贝格积分法：对于复化梯形求积公式，取积分近似值对复化辛普森求积公式因为，所以上述两式相除得所以，同理，对，和分析可得龙贝格积分算法如下：如果 则取，否则，继续计算直到满足条件。程序框图：程序使用说明:运行本程序的时候，直接按照提示输入所求积分的原函数 （比如 ），然后按提示依次输入积分下限 ，积分上限和积分精度 ，然后程序便可计算出原函数在 之间的积分数值。算例计算结果：6.2 4. 四阶龙格-库塔法求解常微分方程的初值问题算法原理：用标准4级4阶R-K法求解一阶常微分方程的算法公式为：求解m阶常微分方程将其转化为一阶微分方程组来求解。为此，引进新的变量 , , ,即可将m阶微分方