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

计算方法A上机大作业
张晓璐硕4011班学号：3114009097

共轭梯度法求解线性方程组
算法原理：由定理3.4.1可知系数矩阵A是对称正定矩阵的线性方程组Ax=b的解与求解二次函数极小点具有等价性，所以可以利用共轭梯度法求解的极小点来达到求解Ax=b的目的。
共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式：

产生的迭代序列在无舍入误差假定下，最多经过n次迭代，就可求得的最小值，也就是方程Ax=b的解。
首先导出最佳步长的计算式。
假设迭代点和搜索方向已经给定，便可以通过的极小化

来求得，根据多元复合函数的求导法则得:

令，得到:
，其中
然后确定搜索方向。给定初始向量后，由于负梯度方向是函数下降最快的方向，故第一次迭代取搜索方向。令

其中。第二次迭代时，从出发的搜索方向不再取，而是选取，使得与是关于矩阵A的共轭向量，由此可求得参数:

然后从出发，沿进行搜索得到
设已经求出，计算。
令，选取，使得和是关于A的共轭向量，可得:

具体编程计算过程如下：
给定初始近似向量以及精度；
计算，取；
Fork=0ton-1do
（i）；
（ii）；
（iii）；
（iv）若，则输出近似解，停止；否则，转（v）；
（v）；
（vi）；
Enddo
程序框图：

程序使用说明：本共轭梯度法求解线性方程的程序是一个函数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)；
程序如下：
clearAbx0%程序运行前首先清除A，b和X0的数值，以免影响计算
clc
n=input('请输入对称正定矩阵A的阶数n=');
epsa=input('请输入误差=');
fork=1:(n-1)%读取题目3.2的矩阵A
A(k,k)=-2;
A(k,k+1)=1;
A(k+1,k)=1;
end
A(n,n)=-2;
A;
b(1,1)=-1;%读取题目3.2的矩阵b
b(n,1)=-1;
b;
x0(1,1)=input('假定初始向量各元素相同，均为：');%给题目3.2的迭代过程赋初值
forkk=2:n
x0(kk,1)=x0(1,1);
end
x=Gongetidu2(A,b,x0,epsa)%调用共轭梯度法求线性方程组的函数

functionx=Gongetidu2(A,b,x0,epsa)
n=size(A,1);
x=x0;
r=b-A*x;
d=r;
fork=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]
三次样条差值
算法原理（三次样条插值函数的导出）：
(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)可得三次样条插值函数的数学表达式。
第二种边界条件：已知。记，则有
所以：

即

其中

所以得到第二种边界条件下的三弯矩方程组：

该方程组系数矩