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

1、求得参数0：计算方法 A 上机大作业1.共轲梯度法求解线性方程组算法原理：由定理 3.4.1 可知系数矩阵 A 是对称正定矩阵的线性方程组 Ax=b1 1的解与求解二次函数 f(x)-xf(x)-xTAxbAxbTx x 极小点具有等价性，所以可以利用1 11共腕梯度法求解 f(x)-xAxbxf(x)-xAxbx 的极小点来达到求解 Ax=b 的目的。共腕梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式：(k1)(k)(k)xxxxkd d产生的迭代序列x(1),x(2),x，在无舍入误差假定下，最多经过 n 次迭代，就可求得 f(x)f(x)的最小值，也就是方程 Ax=b

2、的解。首先导出最佳步长k的计算式。假设迭代点 x x(k)和搜索方向 d d(k)已经给定，便可以通过()f(x(k)d(k)的极小化min()f(x(k)d(k)来求得，根据多元复合函数的求导法则得：1()f(x(k)d(k)Td(k)r(k)T(k)-_d-苴中r(k)bAx(k)d(k)T(k)，八十rbAx然后确定搜索方向 d d(k)。给定初始向量 x x(0)后，由于负梯度方向是函数下降最快的方向，故第一次迭代取搜索方向d(0)r(0)f(x(0)bAx0)。令(1)(0),0 xxxx0d dr(0)Td(0)其中0F 一不。第二次迭代时，从x(1)出发的搜索方向不再取r，而d(

3、0)TAd是选取 d dr rd d(0),使得 d d与 d d是关于矩阵 A 的共腕向量，由此可r(1)TAd0d(0)TAd(0)然后从 x x(1)出发，沿 d d(1)进行搜索得到 X Xx xidid设已经求出 X X(k1)X X 的kd d(k)，计算 r r(k1)bAxbAx(k1)o令 d d(k1)r r(k1)kd d(k)，选取k,使得 d d(k1)和 d d(k)是关于 A 的共腕向量，可得：r(k1)TAd(k)kd(k)TAd(k)具体编程计算过程如下：(1)给定初始近似向量x(0)以及精度0;(2)(2)计算 r rbAxbAx，取 d dr r；(3)F

4、ork=0ton-1dor(k)T(k)rd00d(k)TAd(k)(ii)x(k1)/kd(k);(k1)(k1)(iii)rbAxrbAx；(iv)(iv)若|r r(k1)(k1)| |或女中-1 1,则输出近似解 x x(k1),停止；否则，转(v);(v);Enddo(k1)r(k)r122(k1)*1)(vi)dr(k)kd；程序框图:输出近似解x(k1)(k1)程序使用说明：本共腕梯度法求解线性方程的程序直接打开 matlab 运行,在求解线性方程组 Ax=b(A 是对称正定矩阵)的时候，直接运行程序 Gongetidufa,输入A,b 的值，虽然该函数是通用的，但是对于矩阵 A

5、 和向量 b 的输入，需要使用者根据A 和 b 的特点自行输入。算例 3.4.2 计算结果：对 99 页例题 3.4.2,运行程序 Gongetidufa 将矩阵 A,b 读入系统程序如下：clearallclcA=input(请输入 A 的值)；输入 n 阶正定矩阵刖勺值b=input(请输入 b 的值:);%输入 b 的值n=size(A,1);%求出矩阵 A 的行数x=zeros(n,1);%给定 x 的初始值e=10A(-12);%给出精度r=b-A*x;d=r;for(i=1:n)a=norm(r,2)A2/(d*A*d);%求最佳步长x=x+a*d;j=r；r=b-A*x;if(n

6、orm(r)=e|i=n)break;elseB=norm(r,2)A2/norm(j,2)A2;d=r+B*d;end2,三次样条差值算法原理(三次样条插值函数的导出)：(i) .导出在子区间xi1，x上的 S(x)的表达式由于 S(x)的二阶导数连续，设 S(x)再节点 X 处的二阶导数值 S(xi)=Mi,其中 Mi 为未知的待定参数。由 S(x)是分段的三次多项式知，S(x)是分段线性函数，S(x)在子区间为1?上可表示为S，(x)-MMixi1xiXX1xixxx1一一 MMi1MMi,xi1xxi。hi其中 hi=xi-x(i-1),对上式两次积分得到hi2h2b(yi1-Mi1)

7、/h,Ci(yi+“)加66将b和q代入 S(x)S(x)可得endx计算结果：x=1;1;1%俞出最终的 x 的结果3xixS(x)-M6nbi(Xx)3xxi1G(x61x1),Mixi1xxi由插值条件S(X1)y,S(x)yi得到(ii) .建立参数M的方程组上式中令xx得 S(x)在 xi 处的左导数 S(x)S(x)一.一一_1,._1.因为 S(x)在内节点 xi 处一阶导数连续，所以 S(x)SS(x)S(X),(X),进一步推导可得M12MiM1di,i1,2,3,.n1其中hii,hh1di6(皿，幺上)6fX1,x,Xi1,i1,2,.,n1hih1h1hi上式为三弯矩方

8、程组，因为三弯矩方程组只有 n-1 个方程，不能确定 n+1 个未知量Mi,所以需要再增加两个方程，由边界条件确定。第一种边界条件：此时已知f(a)和f(b).不妨取MOfa),Mnf”b),这时三弯矩方程组化为：2MI1M2d1MiMi12MiiMi1dii2,3,.,n2n1Mi22Mn1dn1nM以上方程组系数矩阵式严格三对角占优矩阵，可用追赶法求解。求出Mi(i1,2,.n1)后，代入 S(x)可得三次样条插值函数的数学表达式。第二种边界条件：已知f(a)和f(b)。记YOf(a),Ynf(b),则有S(x)XiX6KMi1 1XxiM6hi(yi1hi2 23M(XiX)hi2 2(

9、Yi-6Mi)/hi(XXi1),Xi对S(X)求导可得S(X)xxi2hiMiMi16hi,Xi12hi2MiXXi(YiYi1)/hi4 一一一，令xX1得到右导数 S(X),S(X),S(Xo)yo,S(Xn)yn所以：其中所以得到第二种边界条件下的三弯矩方程组:2MoMido,iM2MiMrdi1,2,3,.,n1,Mni2Mndn该方程组系数矩阵是严格三对角占优矩阵，可用追赶法求解，具体追赶法的求解过程见数值分析教材。第三种边界条件：周期型边界条件.已知 yf(x)yf(x)是以Tbaxn为为周期的周期函数，则由周期性可知，yyo,yniyi,MnMo,MniMi,几i几，这时，可以

10、将点Xn看成内点，则方程组对 i=n 也成立，既有nMni2MnnMnidn,也即nMiMi2Mndn,其中,6vynynyn1dn()hnhlhlhn于是三弯矩方程组化为yiy。hi,hnyo,-M6hnn1Mn3ynyn1,hnyn n，2M0M1doMni2Mndndo6,ViV。yo)dn2(yhnynnyn1)I2MMiMna,iMi12MiiMi1d，i2,3,4,.n,1,nMinMniMB该方程组可用 matlab 直接解出。程序框图如下：开始输入x,Vnsiz0 x2)h(k)x(k)x(k1),k2,3n(k)h(k)/(h(k)h(k1),k2,3,.n11d(k)6(y

11、(k1)y(k)/h(k1)(V(k)V(k1)yh(k)/(h(k)h(k1),k2,3,.n1输入边界条件输出(Xx)3(xX1)3hXxS(x)%JMi1i1.Mi(V1:Mj,56h6h程序使用说明：本程序是求解 137 页例题 4.6.1 的运行结果，通过程序便可求得M,然后根据S(x)3M(2L2Mj(yyMr)J6hi6hi6h2/hixXii(yiMi),Xi1xXi6hi便可得到，在XiXXi上的三次样条插值函数S(x),进而得到整个区间上的三次样条差值函数 S(x)S(x)。算例计算结果：137 页例题 4.6.1 的计算实习1、打开 matlab 运行 Sanciyang

12、tiao 程序2、自行输入 x 和 y 的节点值3、得出计算结果清输入K的值二-2,-L033n4j清潜入描tS：7.11,26,弧293=x3+9+x2+十203*x2-K3+19*x+263*x2-2*xft3+19*x+26-S0*x2+20S*x-1633.龙贝格积分法对于复化梯形求积公式，取积分近似值(TT)T(TT)Tn.(T2nTn) )T2n4141对复化辛普森求积公式baba4h h4f f()()，a a28802880I(f)I(f)I(fI(f) )S Sn一一bahbah41S2n频(万,因为 f f(4)( () )(4),、f f( (1),),所以上述两式相除得

13、所以,I(f)I(f)&)S S2n同理,I(f)Cn) Cn对 G,SG,S2n和 C C2n分析可得sn711(T2n41玉)C2n142113.41(S2nS)GnG)龙贝格积分算法如下:baT-f(a)f(b),22KJ1T2Kbraf(a(2i1)号),k0,1,2.,22i12S2kC2k0LrT2k141(T2k1T2k)，-1_S?k1-2(S2k1Sk),41-1_-、Ck1-o(Ck1Ck),243122,k0,1,2.,如果%1R2k，则取 IfRIfR2k1,否则，继续计算直到满足条件。程序框图:程序使用说明：运行本程序的时候，直接按照提示输入所求积分的原函数f

14、(X)(比如工)，然后按提示依次输入积分下限a,积分上限 b 和积分1X1X精度epS1,然后程序便可计算出原函数f(X)在a，b之间的积分数值。算例计算结果：209 页题 6.2 第一小题计算实习结果11 1dx0.6931dx0.693101 1x x4.四阶龙格-库塔法求解常微分方程的初值问题算法原理：用标准 4 级 4 阶 R-K 法求解一阶常微分方程的算法公式为:1yi1yi-(Ki2K22K3K4)6c 开始)*输入f,abepsl0,1,21S2kT2k1411(T2k1Tk), k0,1,2C2kS2k1421(S2k1S2k),k0,1R2kC2k11431(C2k1C2k),k0baTi2f(a)f(b)1baba2T2k2kif(a1)2ki),kk3Kihf(x,y)11K2hf(Xi-h,y2K1)r11K3hf(x/y5K2)K4hf(Xih,yiK3)程序框图：开始定义变量 x,y;输入 f(x,y)的函数表达式输入积分区间、步长 h、初始值 y。、计