数值分析上机题3.docx
《数值分析上机题3.docx》由会员分享，可在线阅读，更多相关《数值分析上机题3.docx（14页珍藏版）》请在冰豆网上搜索。
数值分析上机题3
数值分析上机题目3
实验一
1．根据Matlab语言特点，描述Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法。
2．编写Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法的M文件。
3．给定
为五对角矩阵
（1）选取不同的初始向量
及右端面项向量b，给定迭代误差要求，分别用编写的Jacobi迭代法和Gauss-Seidel迭代法程序求解，观察得到的序列是否收敛？
若收敛，通过迭代次数分析计算结果并得出你的结论。
（2）用编写的SOR迭代法程序，对于
（1）所选取的初始向量
及右端面项向量b进行求解，松驰系数ω取1<ω<2的不同值，在
时停止迭代，通过迭代次数分析计算结果并得出你的结论。
实验1
1、根据MATLAB语言特点，描述Jacobi迭代法，Gauss-Seidel迭代法和SOR迭代法。
2、编写Jacobi迭代法，Gauss-Seidel迭代法和SOR迭代法的M文件。
Jacobi迭代法
function[x1,k]=GS_2（A,b）
n=length（A）;
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
x1=zeros（n,1）;x0=3*ones（n,1）;k=0;
whilenorm（x1-x0,1）>10^（-7）&k<100
k=k+1;
x0=x1;
x1=D\（（L+U）*x0+b）;
end
k=k
x=x1
Gauss-Seidel迭代法
function[x1,k]=GS_h（A,b）
n=length（A）;
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
x1=zeros（n,1）;x0=3*ones（n,1）;k=0;
whilenorm（x1-x0,1）>10^（-7）&k<100
k=k+1;
x0=x1;
x1=（D-L）\U*x0-D\b;
end
k=k
x=x1
SOR迭代法
function[x1,k]=SOR_h（A,b）
n=length（A）;
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
x1=zeros（n,1）;x0=3*ones（n,1）;k=0;
w=0.96;
whilenorm（x1-x0,1）>10^（-7）&k<100
k=k+1;
x0=x1;
x1=（D-w*U）\（（（1-w）*D+w*L）*x0+w*b）;
end
k=k
x=x1
3、采用Jacobi迭代法，Gauss-Seidel迭代法求解五对角矩阵
clear,clc
A=diag（3*ones（20,1））+diag（（-0.5）*ones（19,1）,-1）+diag（（-0.5）*ones（19,1）,1）+diag（（-0.25）*ones（18,1）,-2）+diag（（-0.25）*ones（18,1）,2）;
b=sum（A'）';
[x1,k1]=Jacob_h（A,b）
[x2,k2]=GS_h（A,b）
运行结果：
两种方法都收敛，k1=27,k2=13。
说明Gauss-Seidel迭代速度比Jacobi迭代速率快
4、采用SOR迭代法程序对五对角矩阵进行求解
clear,clc
A=diag（3*ones（20,1））+diag（（-0.5）*ones（19,1）,-1）+diag（（-0.5）*ones（19,1）,1）+diag（（-0.25）*ones（18,1）,-2）+diag（（-0.25）*ones（18,1）,2）;
b=sum（A'）';
[x3,k3]=SOR_h（A,b）
运行结果当w=0.5时，k3=53，当w=0.96时，k3=19，当w=1.1时，k3=14，当w=1.5时，k3=33。
该结果说明，当w选取合适时，可以大大加快运算速率。
实验二
题目：
多项式最小二乘法
摘要：
对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观察和测量给出离散的一些点，再来求出具体的函数解析式。
又因为测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。
最小二乘法是求解这一问题的很好的方法，本实验运用这一方法实现对给定数据的拟合。
数学原理：
对于给定的测量数据（xi,fi）（i=1,2,…，n），设函数分布为
特别的，取
为多项式
（j=0,1,…，m）
则根据最小二乘法原理，可以构造泛函
令
（k=0,1,…，m）
则可以得到法方程
求该解方程组，则可以得到解
，因此可得到数据的最小二乘解
程序设计：
编写求解多项式拟合的Matlab函数子程序
实验要求：
用最小二乘法处理下面的实验数据.
xi
3
4
5
6
7
8
9
fi
2.01
2.98
3.50
5.02
5.47
6.02
7.05
并作出
的近似分布图。
分别采用一次，二次、五次和偶数次多项式来拟合数据得到相应的拟合多项式，并分别作出它们的曲线图。
实验2
clc
clear
n=input（'请输入多项式拟合次数：
'）;
symst
fori=1:
n
f（i）=t^（i）;
end
x=[3456789]';
y=[2.012.983.505.025.476.027.05]';
m=length（x）;
fori=1:
m
forj=1:
n
hs=inline（f（j）,'t'）;
A（i,j）=hs（x（i））;
end
end
h=ones（m,1）;
A=[hA];
A1=A'*A;
y1=A'*y;
x1=A1\y1;
f=0;
fori=1:
n+1
f=f+x1（i）*t^（i-1）;
end
plot（x,y,'*'）;
holdon
x1=flipud（x1）;
x2=linspace（min（x）,max（x））;
y2=polyval（x1,x2）;
tt=poly2str（x1,'x'）
text（5,7,0,tt）
plot（x2,y2）
实验三
实验名称：
非线性方程组数值求解的Newton类方法试验。
实验目的：
用Newton类方法求解线性方程组F（x）=0，理解其解的复杂性、初始点选择策略、减少算法工作量的方法等。
实验容与要求：
分别用Newton法用Broyden秩1校正法求解下面非线性方程组
（1）写出MATLAB源代码；
（2）给出迭代五次以上的结果；
（3）尝试不同的初值，如可取
）；
（4）计算两种方法的用时。
实验3
1、2、
采用Newton法·¨
clear,clc
x0=[0,0,0]';
y0=f（x0）;
yy0=df（x0）;
x1=x0-yy0\y0;
k=1;formatlong
whilenorm（x1-x0,1）>10^（-5）&k<100
k=k+1;
x0=x1;
y0=f（x0）;
yy0=df（x0）;
x1=x0-yy0\y0;
end
x=x1
k=k
运行结果：
x=
0.0332
-0.5047
-0.2383
k=
45
Broyden秩1法
clear,clc
x0=[0,0,0]';
y0=f（x0）;
A0=df（x0）;
x1=x0-A0\y0;
y1=f（x1）;
k=1;formatlong
whilenorm（x1-x0,1）>10^（-2）&k<100000
k=k+1;
g=y1-y0;
y=x1-x0;
A1=A0+（g-A0*y）/（y'*y）*y';
x0=x1;
x1=x0-A1\y1;
A0=A1;
y0=f（x0）;
y1=f（x1）;
end
x=x1
k=k
运行结果，
x=
0.3973
-0.3592
-0.1240
k=
5
3.尝试不同初值。
x0=[0.1,0.1,0.1]
采用Newton法·¨
x=
0.2459
-0.6348
-0.0848
k=
68
Broyden秩1法
x=
0.0213
-0.4117
-0.2104
k=
7
该结果说明，采用不同的初值，对结果没有影响。
4、计算两种方法的用时
采用Newton法·¨
t1=
0.7623
Broyden秩1法
t2=
0.2444
采用Broyden秩1法运算比较快。
升级会员 