上QQ阅读APP看书，第一时间看更新
1.7 线性方程组求解
线性方程组由n个方程组合而成，一般情况下，方程的个数和未知数的个数相同，这样的线性方程组一般具有唯一解。如果方程的个数小于未知数的个数，大多数情况下，该线性方程组含有较多的可行解。线性方程组求解方法较多，不同的求解方法求解矩阵精度是有差异的，本节着重讲述应用较多的牛顿Bewton迭代法、高斯-赛德尔Gauss-Seidel迭代法以及雅克比Jacobi迭代法。
1.7.1 牛顿Newton迭代法
1．Newton插值
Newton均差插值公式为：
其中，
是k阶均差，可由均差表1-12方便计算得到。
表1-12 均差计算表格
在MATLAB中实现利用均差的牛顿插值的代码如下：
function f = Newton_interp(x,y,x0)
% x,y为输入输出变量
% x0为插值点处的横坐标值
syms t; % 符号变量
if(length(x) == length(y)) % 维数相等，说明数据维数合理
n = length(x); % 数据长度
c(1:n) = 0.0; % 初始化
else
disp('x和y的维数不相等！');
return;
end
f = y(1); % 初始值
y1 = 0; % 初始化
l = 1; % 初始化
for i = 1:(n-1)
for j = (i+1):n
y1(j) = (y(j)-y(i))/(x(j)-x(i)); % 斜率
end
c(i) = y1(i+1); % 斜率值
l = l*(t-x(i)); % 均差
f = f + c(i)*l; % 直线方程
simplify(f); % 简化f直线方程
if(i==n-1)
if(nargin == 3)
f = subs(f,'t',x0); % 赋值，在x0处的y值
else
f = collect(f); % 将插值多项式展开
f = vpa(f, 6); % 显示精度为6位小数
end
end
y = y1; % 赋值
end
如
上各阶导数存在，则在此区间取n个节点进行Newton插值，主函数编程如下：
% Designed by Yu Shengwei From SWJTU University
% 2014年12月29日
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
x = -5:1:5; % x
y = 1./(1+x.*x); % y
% 精确解曲线
xj = -5:0.1:5; % 待插值的x
yj = 1./(1+xj.*xj); % 精确解
plot(xj,yj,'linewidth',2) % 画图
hold on % 图形保持句柄
% 高次多项式插值
for i = 1:length(xj)
disp([strcat('当前插值点数为 ',num2str(i)),strcat('总插值点数为 ',num2str(length(xj)))])
yc(i) = Newton_interp(x,y,xj(i));
end
plot(xj,yc,'r--','linewidth',2)
grid on;xlabel('x'),ylabel('y')
legend('原数据曲线','牛顿插值曲线')
toc % 计时结束'
运行程序输出结果如图1-42所示。
图1-42 Newton插值
2．Newton迭代法
将非线性方程线性化，以线性方程的解逐步逼近非线性方程的解，这就是Newton迭代法的基本思想。
设已知方程f（x）＝0的近似根x0，f（x）在其零点x*邻近一阶连续可微，且
，当x0充分接近x*时，f（x）可用Taylor公式近似表示为：
取此x作为原方程的新近似值x1，重复以上步骤，于是得迭代公式：
在MATLAB中编程实现的牛顿法的函数程序代码如下：
function xn=Newton_dd(fun,x0,D)
% 牛顿法求解非线性方程的解
% 函数输出
% xn：为所求非线性方程的解
% 函数输入
% fun：为所定义的函数
% x0：为初始值
% D：为计算的精确度
[f0,df]=feval(fun,x0); % 一阶导数
if df==0;
error('d[f(x)/dx]=0 at x0'); % 一阶导数为0，不满足牛顿迭代
end
if nargin<3;
D=1e-6; % 默认值
end
d=f0/df; % 初始斜率
while abs(d)>D; % 判断条件
x1=x0-d; % 解的更新
x0=x1; % 赋值
[f0,df]=feval(fun,x0); % 函数值和一阶导数
if df==0; % 一阶导数
error('d[f(x)]/dx=0 at x0'); % 一阶导数为0，不满足牛顿迭代
end
d=f0/df; % 更新斜率
end
xn=x1; % 赋值
end
利用牛顿法求方程x-lnx＝2于区间［2，4］的根。
首先定义计算函数的函数文件如下：
function [y,dy,d2y] = Newton_fun_diff(x)
% y为在x处的函数值
% dy为一阶导数
% d2y为二阶导数
y = x-log(x)-2; % 计算函数值
if nargout>1;
ff=sym('x-log(x)-2'); % 定义符号函数
dy=diff(ff); % 求一阶导数
dy=subs(dy,x); % 赋值
end
if nargout==3;
d2y=diff(ff,2); % 求二阶导数
d2y=subs(d2y,x); % 赋值
end
由此编写牛顿迭代算法编程如下：
% Designed by Yu Shengwei From SWJTU University
% 2014年12月29日
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
x = linspace(2,4,200); % 对自变量取样
y = Newton_fun_diff(x); % 计算各点的函数值
plot(x,y,'linewidth',2); % 绘制函数y(x)曲线
hold on;
plot(xlim,[0,0],'m:'); % 绘制零刻度线
x0 = 3; % 初始化解的位置
xj=Newton_dd('Newton_fun_diff',x0) % 牛顿求根，初始点是3
plot(xj,Newton_fun_diff(xj),'rs'); % 绘制解对应的点
toc % 计时结束'
运行程序输出结果如下：
xj =
(14176964281861263*log(14176964281861263/4503599627370496))/9673364654490767 + 14176964281861263/9673364654490767
ans =
3.14619
时间已过 0.325242 秒。
输出图形如图1-43所示。
图1-43 牛顿迭代法求解结果
1.7.2 高斯-赛德尔Gauss-Seidel迭代法
高斯-塞德尔（Gauss-Seidel）迭代格式如下：
其分量形式为：
采用高斯-赛德尔迭代法进行求解，将系数矩阵A分解为A＝D–L–U，其中：
其矩阵形式如下：
在MATLAB中编程实现Gauss-Seidel迭代法一般化程序如下：
function x = gauss_seidel_x(A,B,x0,Err)
% A为方程组系数
% B为方程组值
% x0为初值
% Err为求解精度
D = diag(diag(A)); % 提取A中对角元素
L = -tril(A)+D; % 求下三角矩阵
U = -triu(A)+D; % 求上三角矩阵
DL = D-L;
A_DL = inv(DL); % 求逆
x = A_DL*U*x0+A_DL*B;
while norm(x-x0)>Err % d当两次计算结果2范数小于Err退出循环体
x=x0;
x0 = A_DL*U*x+A_DL*B;
end
用Gauss-Seidel迭代法求解下列方程组：
由上述方程组可得
，设初值
。
采用Gauss-Seidel迭代法进行求解，程序如下：
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
format shortE
tic % 运算计时
A = [6,2,-1;
1,4,-2;
-3,1,4]; % 方程系数矩阵
b = [-3,2,4]';
x0 = [-0.5,0.5,0.25]';
% gauss_seidel求解
y = gauss_seidel_x(A,b,x0,eps)
toc % 计时结束
运行程序输出结果如下：
y =
-7.2727e-01
8.0808e-01
2.5253e-01
时间已过 0.571006 秒。
1.7.3 雅克比Jacobi迭代法
设有n阶线性方程组：
由此可建立迭代格式：
简记为：
在MATLAB中编程实现Jacobi迭代法的程序代码如下：
function [x,n]=jacobi(A,b,x0,eps,varargin)
% 采用Jacobi迭代法求线性方程组Ax=b的解
% 函数输入
% A：线性方程组的系数矩阵
% b：线性方程组中的常数向量
% x0：迭代初始向量
% eps：解的精度控制
% varargin：迭代步数控制
% 函数输出
% x：线性方程组的解
% n：求出所需精度的解实际的迭代步数
if nargin==3
eps= 1.0e-6; % 默认值
elseif nargin<3
error('输入错误') % 输入错误
return
elseif nargin ==5
M = varargin{1};
end
M = 500; % 迭代次数
D=diag(diag(A)); % 求A的对角矩阵
L=-tril(A,-1); % 求A的下三角矩阵
U=-triu(A,1); % 求A的上三角矩阵
B=D\(L+U); % 矩阵运算
f=D\b; % 矩阵相除
x=B*x0+f;
n=1; % 迭代次数
% 迭代过程
while norm(x-x0)>=eps
x0=x; % 赋值
x =B*x0+f; % 值更新
n=n+1;
if(n>=M)
disp('Warning:迭代次数太多，可能不收敛！');
return;
end
end
对于下列方程组：
采用雅克比矩阵迭代法求解，编程如下：
% Designed by Yu Shengwei From SWJTU University
% 2014年12月29日
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
format shortE
tic % 运算计时
A=[9,-1,-1; % 系数矩阵
-1,8,0;
-1,0,9];
b=[7,7,8]'; % 方程组右边系数
x0 = [0,0.1,0.1]';
[x,n]=jacobi(A,b,x0,eps) % 采用Jacobi迭代法求线性方程组Ax=b的解
toc % 计时结束
运行程序输出结果如下：
x = % 方程的解
1
1
1
% 迭代次数
n =
22
时间已过 0.002594 秒。