迭代解法
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
迭代解法
迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中,迭代解法主要包括Jacobi 迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。
1.Jacobi迭代法
对于线性方程组Ax=b,如果A为非奇异方阵,即aii≠0(i=1,2,…,n),则可将A分解为A=D-L-U,其中D为对角阵,其元素为A的对角元素,L与U为A的下三角阵和上三角阵,于是Ax=b化为:
x=D\(L+U)x+D\b
与之对应的迭代公式为:
x(k+1)=D\(L+U)x(k)+D\b
这就是Jacobi迭代公式。如果序列{x(k+1)}收敛于x,则x必是方程Ax=b的解。
Jacobi迭代法的MATLAB函数文件Jacobi.m如下:
function [y,n]=jacobi(A,b,x0,eps)
if nargin==3
eps=1.0e-6;
elseif nargin<3
error
return
end
D=diag(diag(A)); %求A的对角矩阵
L=-tril(A,-1); %求A的下三角阵
U=-triu(A,1); %求A的上三角阵
B=D\(L+U);
f=D\b;
y=B*x0+f;
n=1; %迭代次数
while norm(y-x0)>=eps
x0=y;
y=B*x0+f;
n=n+1;
end
例7-5 用Jacobi迭代法求解下列线性方程组。设迭代初值为0,迭代精度为10-6。
在命令中调用函数文件Jacobi.m,命令如下:
A=[10,-1,0;-1,10,-2;0,-2,10];
b=[9,7,6]';
[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)
2.Gauss-Serdel迭代法
在Jacobi迭代过程中,计算时,已经得到,不必再用,即原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b,于是得到:
x(k+1)=(D-L)\Ux(k)+(D-L)\b
该式即为Gauss-Serdel迭代公式。和Jacobi迭代相比,Gauss-Serdel迭代用新分量代替旧
分量,精度会高些。
Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下:
function [y,n]=gauseidel(A,b,x0,eps)
if nargin==3
eps=1.0e-6;
elseif nargin<3
error
return
end
D=diag(diag(A)); %求A的对角矩阵
L=-tril(A,-1); %求A的下三角阵
U=-triu(A,1); %求A的上三角阵
G=(D-L)\U;
f=(D-L)\b;
y=G*x0+f;
n=1; %迭代次数
while norm(y-x0)>=eps
x0=y;
y=G*x0+f;
n=n+1;
end
例7-6 用Gauss-Serdel迭代法求解下列线性方程组。设迭代初值为0,迭代精度为10-6。在命令中调用函数文件gauseidel.m,命令如下:
A=[10,-1,0;-1,10,-2;0,-2,10];
b=[9,7,6]';
[x,n]=gauseidel(A,b,[0,0,0]',1.0e-6)
例7-7 分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性方程组,看是否收敛。
命令如下:
a=[1,2,-2;1,1,1;2,2,1];
b=[9;7;6];
[x,n]=jacobi(a,b,[0;0;0])
[x,n]=gauseidel(a,b,[0;0;0])
四阶龙格库塔法
M文件如下
%xout tout 为输出
%func为要求的函数方程
%ts为求值范围和精度
%x0为初始值
下面的是名为rk_4的主程序
function[tout,xout]=rk_4(func,ts,x0) t0=ts(1);
tf=ts(2);
if length(ts)==3
h=ts(3);
else
h=(tf-t0)/100;
end
xout=x0;
tout=t0;
for t=t0:h:tf-h
k1=func(t,x0);
k2=func(t+h/2,x0+k1*h/2);
k3=func(t+h/2,x0+k2*h/2);
k4=func(t+h,x0+k3*h);
x0=x0+(k1+2*k2+2*k3+k4)*h/6;
xout=[xout;x0];
tout=[tout;t+h];
end
要求的函数要写成函数文件如下
function dy=newfunction(t,y)
dy(1)=4*y(1)-2*y(1)*y(2);
dy(2)=y(1)*y(2)-3*y(2);
在命令窗口写的程序为
[xout,tout]=rk_4(@newfunction,[0,9,0.1],[2,3])
输出为xout,tout的矩阵,若要画图
>> plot(xout,tout)
这就是一道题目的完整解,要解的方程为