上机作业.docx
《上机作业.docx》由会员分享，可在线阅读，更多相关《上机作业.docx（39页珍藏版）》请在冰豆网上搜索。
上机作业
上机作业12
1、分别用不动点迭代与Newton法求解方程3x-5^x+4=0的正根与负根2
2、用Newton法与重根计算法求解方程x-sinx＝0的根.再用Steffensen’smethod加速Newton法收敛，比较结果5
上机作业28
分别用Jacobi、Seidel、Sor（w=0.8,1.1,1.2,1.3,1.4,1.5）方法求解方程组Ax=b8
上机作业320
用Newton法与最速下降法求方程组x2+2y2-2=0,x2=y在（0.8，0.7）附近的根20
上机作业425
1、用幂法与反幂法求矩阵A的按模最大、最小特征值与对应的特征向量25
2、用Householder变换求矩阵A的QR分解，并用QR方法做3次迭代25
上机作业527
目的：
观察lagrange插值的Runge现象27
上机作业627
目的：
掌握数值求积公式27
上机作业1
1、分别用不动点迭代与Newton法求解方程3x-5^x+4=0的正根与负根。
解：
1）实验原理
不动点迭代：
（1）构造不动点方程x=ϕ（x），即由方程f（x）=0等价变形为x=ϕ（x）（这样，求f（x）的零点可转化成求ϕ（x）的不动点），ϕ（x）称为迭代函数;
（2）建立迭代格式，xk+1=ϕ（xk）称为不动点迭代法。
当给定初值x0后,由迭代格式可求得数列{xk}。
如果{xk}收敛于x*，且ϕ（x）在x*连续，则x*就是方程的根。
2）不动点迭代方程：
由方程3x-5^x+4=0，构造迭代方程x=（5^x-4）/3求负根，x=log（3x+4）/log5求正根。
给定初值：
p0=1.5,p1=-1.5.
Newton迭代：
设r是f（x）=0的根，选取x0作为r初始近似值，过点（x0,f（x0））做曲线y=f（x）的切线L，L的方程为y=f（x0）+f'（x0）（x-x0），求出L与x轴交点的横坐标x1=x0-f（x0）/f'（x0），称x1为r的一次近似值。
过点（x1,f（x1））做曲线y=f（x）的切线，并求该切线与x轴交点的横坐标x2=x1-f（x1）/f'（x1），称x2为r的二次近似值。
重复以上过程，得r的近似值序列，其中x（n+1）=x（n）－f（x（n））/f'（x（n）），称为r的n+1次近似值，上式称为牛顿迭代公式。
Newton迭代方程：
x=x-（3*x-5^x+4）/（3-（5^x）*ln5）;给定初值：
p0=1.5,p1=-1.5.
3）程序及结果
不动点迭代程序：
disp（'不动点迭代法'）;
n0=100;p0=-1.5;
fori=1:
n0
p=（5^p0-4）/3;
ifabs（p-p0）<=10^（-6）
ifp<0
disp（'|p-p0|='）
disp（abs（p-p0））
disp（'不动点迭代法求得方程的负根为:
'）
disp（p）;
disp（'用不动点迭代法迭代步数为：
'）
disp（i）;
break;
else
p0=p;
end
end
p1=1.5;
fori=1:
n0
p=（log（3*p1+4））/log（5）;
ifabs（p-p1）<=10^（-6）
ifp>0
disp（'|p-p1|='）
disp（abs（p-p1））
disp（'用不动点迭代法求得方程的正根为：
'）
disp（p）;
disp（'用不动点迭代法迭代步数为：
'）
disp（i）;
break;
else
p1=p;
end
end
Newton迭代程序：
disp（'牛顿法'）
n0=80;p0=1.5;
fori=1:
n0
p=p0-（3*p0-5^p0+4）/（3-log（5）*（5^p0））;
ifabs（p-p0）<=10^（-6）
ifp>0
disp（'|p-p0|='）
disp（abs（p-p0））
disp（'用牛顿法求得方程的正根为:
'）
disp（p）;
disp（'用不动点迭代法迭代步数为：
'）
disp（i）;
end
break;
else
p0=p;
end
end
p1=-1.5;
fori=1:
n0
p=p1-（3*p1-5^p1+4）/（3-log（5）*（5^p1））;
ifabs（p-p1）<=10^（-6）
ifp<0
disp（'|p-p1|='）
disp（abs（p-p1））
disp（'用牛顿法求得方程的负根为:
'）
disp（p）;
disp（'用不动点迭代法迭代步数为：
'）
disp（i）;
end
break;
else
p1=p;
end
end
ifi==n0
disp（n0）
disp（'用牛顿迭代后无法求出方程的解'）
end
运行结果：
结果分析：
选用不动点迭代公式x=（5^x-4）/3求负根的运行结果是收敛的，求正根发散，这是由于指数函数的存在使迭代发散。
而使用的迭代公式x=log（3x+4）/log5在求负根时发散，求正根时收敛，这是由于公式中存在对数函数。
由以上运算结果可知，Newton迭代法公式唯一，且求正负根均收敛，迭代步数明显少于不动点迭代，运算效率高。
2、用Newton法与重根计算法求解方程x-sinx＝0的根.再用Steffensen’smethod加速Newton法收敛，比较结果。
解：
Newton迭代法原理：
设r是f（x）=0的根，选取x0作为r初始近似值，过点（x0,f（x0））做曲线y=f（x）的切线L，L的方程为y=f（x0）+f'（x0）（x-x0），求出L与x轴交点的横坐标x1=x0-f（x0）/f'（x0），称x1为r的一次近似值。
过点（x1,f（x1））做曲线y=f（x）的切线，并求该切线与x轴交点的横坐标x2=x1-f（x1）/f'（x1），称x2为r的二次近似值。
重复以上过程，得r的近似值序列，其中x（n+1）=x（n）－f（x（n））/f'（x（n）），称为r的n+1次近似值，上式称为牛顿迭代公式。
重根计算法原理：
在Newton迭代法的基础上改进迭代公式，使得x（n+1）=x（n）－mf（x（n））/f'（x（n）），m为根的重数。
牛顿法：
Newton迭代方程：
x=x-（x-sinx）/（1-cosx）,给定初值：
x=0
程序：
disp（'牛顿法'）
n0=80;p0=1;
fori=1:
n0
p=p0-（p0-sin（p0））/（1-cos（p0））;
ifabs（p-p0）<=10^（-6）
disp（'|p-p0|='）
disp（abs（p-p0））
disp（'用牛顿法求得方程的根为:
'）
disp（p）;
disp（'用Newton法迭代步数为：
'）
disp（i）;
break;
else
p0=p;
end
end
ifi==n0
disp（n0）
disp（'用牛顿迭代后无法求出方程的解'）
end
重根计算法程序：
disp（'重根计算法'）
n0=80;p0=1;m=2;
fori=1:
n0
p=p0-m*（p0-sin（p0））/（1-cos（p0））;
ifabs（p-p0）<=10^（-6）
disp（'|p-p0|='）
disp（abs（p-p0））
disp（'用重根计算法求得方程的根为:
'）
disp（p）;
disp（'用重根计算法迭代步数为：
'）
disp（i）;
break;
else
p0=p;
end
end
Steffensen加速Newton法程序：
disp（'Steffensen加速牛顿法'）
n0=80;p0=1;i=0;p
（1）=p0;
whilei<=n0
fork=1:
2
p（k+1）=p（k）-（p（k）-sin（p（k）））/（1-cos（p（k）））;
end
p1=p
（1）-（p
（2）-p
（1））^2/（p（3）-2*p
（2）+p
（1））;
p=p1-（p1-sin（p1））/（1-cos（p1））;
ifabs（p-p1）<=10^（-6）
disp（'|p-p1|='）
disp（abs（p-p1））
disp（'用Steffensen加速牛顿法求得方程的根为:
'）
disp（p）;
disp（'用Steffensen加速Newton法迭代步数为：
'）
i=i+1;
disp（i+1）;
break;
else
p
（1）=p1;
end
end
运行结果：
结果分析：
由运行结果可知，Newton法收敛较慢，重根计算法效率大大提高，而Steffensen加速牛顿法效率最高，仅两次迭代就达到了要求的计算精度。
上机作业2：
分别用Jacobi、Seidel、Sor（w=0.8,1.1,1.2,1.3,1.4,1.5）方法求解方程组Ax=b。
这里
X0=[0……0]’;e=10-5
要求：
抄题，公式，程序、计算结果（终止迭代步数k、近似解x（k）），结果分析（三种迭代列表）。
实验原理：
1、雅可比迭代法（J-迭代法）：
线性方程组
，可以转变为：
迭代公式
其中
，称
为求解
的雅可比迭代法的迭代矩阵。
以下给出雅可比迭代的分量计算公式，令
，由雅可比迭代公式有
，既有
，于是，解
的雅可比迭代法的计算公式为
2、高斯-赛德尔迭代法（GS-迭代法）：
GS-迭代法可以看作是雅可比迭代法的一种改进，给出了迭代公式：
其余部分与雅克比迭代类似。
3、逐次超松弛迭代法（SOR-迭代法）：
选取矩阵A的下三角矩阵分量并赋予参数w，将之作为分裂矩阵M，
，其中，w>0，为可选择的松弛因子，又
（1）公式构造一个迭代法，其迭代矩阵为
从而得到解
的逐次超松弛迭代法。
其中：
由此，解
的SOR-迭代法的计算公式为
（2）
观察
（2）式，可得结论：
当w=1时，SOR-迭代法为J-迭代法。
当w>1时，称为超松弛迭代法，当w<1时，称为低松弛迭代法。
Jacobi迭代法程序：
disp（'用Jacobi迭代法求解线性方程组'）;
i=1;
i<=100;
a=[-13111111111;1-1311111111;11-131111111;111-13111111;1111-1311111;11111-131111;111111-13111;1111111-1311;11111111-131;111111111-13];
d=diag（diag（a））;
l=d-tril（a）;
u=d-triu（a）;
d0=inv（d）;
b=[-1;-1;-1;-1;-1;-1;-1;-1;-1;-1];
x0=zeros（10,1）;
B=-d0*（l+u）;
f=d0*b;
x=B*x0+f;
whilenorm（x-x0,inf）>=1e-5
x0=x;x=B*x0+f;i=i+1;
end
disp（'方程组的解为:
'）;
x
disp（'迭代次数:
'）
i
break;
Gauss-Seidel迭代法程序：
disp（'用Gauss-Seidel迭代法求解线性方程组'）;
i=1;
i<=100;
a=[-13111111111;1-1311111111;11-131111111;111-13111111;1111-1311111;11111-131111;111111-13111;1111111-1311;11111111-131;111111111-13];
d=diag（diag（a））;
l=d-tril（a）;
u=d-triu（a）;
b=[-1;-1;-1;-1;-1;-1;-1;-1;-1;-1];
x0=zeros（10,1）;
B=-inv（d+l）*u;
f=inv（d+l）*b;
x=B*x0+f;
whilenorm（x-x0,inf）>=1e-5
x0=x;x=B*x0+f;i=i+1;
end
disp（'方程组的解为:
'）;
x
disp（'迭代次数:
'）
i
break;
Sor迭代法程序：
分别令w=0.8,1.1,1.2,1.3,1.4,1.5
disp（'用Sor迭代法求解线性方程组'）;
w=0.8;
i=1;
i<=100;
a=[-13111111111;1-1311111111;11-131111111;111-13111111;1111-1311111;11111-131111;111111-13111;1111111-1311;11111111-131;111111111-13];
d=diag（diag（a））;
l=d-tril（a）;
u=d-triu（a）;
b=[-1;-1;-1;-1;-1;-1;-1;-1;-1;-1];
x0=zeros（10,1）;
B=inv（d+w*l）*[（1-w）*d-w*u];
f=w*inv（d+w*l）*b;
x=B*x0+f;
whilenorm（x-x0,inf）>=1e-5
x0=x;
x=B*x0+f;
i=i+1;
end
disp（'松弛因子为:
'）;
w
disp（'方程组的解为:
'）;
x
disp（'迭代次数:
'）;
i
break;
运行结果：
运行结果分析：
由运行结果截屏的对比可发现，Jacobi迭代法收敛速度较慢，seidel
迭代法与松弛因子为1.1的Sor法的迭代效率一样，而w=0.8时，Sor的效率最高，随着松弛因子的增大，计算效率下降。
上机作业3
用Newton法与最速下降法求方程组x2+2y2-2=0,x2=y在（0.8，0.7）附近的根。
要求：
1）抄题；2）迭代公式（初值）或简单原
理；3）程序，结果；（打印）4）结果分析。
牛顿迭代法原理
对于非线性方程
在x（k）处按照多元函数的泰勒展开，并取线性项得到
其中
初值取x=0.8,y=0.7;
实验步骤
步骤一：
编写牛顿迭代法的基本程序。
打开Editor编辑器，输入以下语句：
function[x,n,data]=new_ton（x0）
ifnargin==1
tol=1e-6;
end
x1=x0-f1（x0）/df1（x0）;
n=1;
while（norm（x1-x0）>tol）
x0=x1;
x1=x0-f1（x0）/df1（x0）;
n=n+1;
data（:
n）=x1;
end
x=x1;
以文件名new_ton.m保存。
步骤二：
编写方程函数与方程的Jacobi矩阵函数。
（1）打开Editor编辑器输入以下语句：
%牛顿迭代法的方程函数
functionf=f1（x0）
x=x0
（1）;
y=x0
（2）;
f1=x^2-2*x-y+0.5;
f2=x^2+4*y^2-4;
%最后方程函数以行向量输出
f=[f1f2];
以文件名f1.m保存。
（2）新打开Editor编辑器输入以下语句：
%牛顿迭代法的jacobi矩阵
functionf=df1（x0）;
x=x0
（1）;
y=x0
（2）;
f=[2*x4*y
2*x-1];
以文件名df1.m保存。
步骤三：
编写主函数。
打开Editor编辑器输入以下语句：
%牛顿迭代法的主函数
x0=[0.80.7];
[x,n,data]=new_ton（x0）;
disp（'计算结果为'）
x
disp（'迭代次数为'）
n
以文件名new_main.m保存。
最速下降法原理
对于非线性方程组
令
如果给定一个初值
，我们希望找到一条路线每一次
迭代以后代价函数都会比原来小一些。
l称为步长因子
，
的不同，就构成了不同的下降算法。
如果取
就是所谓的最速下降法。
最速下降法是大范围收敛的h在某
出沿最速下降方向
下降的最快
实验步骤
步骤一：
打开Editor编辑器，输入以下语句：
symsxy
f1=x^2+2*y^2-2;
f2=x^2-y;
h=f1^2+f2^2
grad=[diff（h,x）,diff（h,y）];
grad=simple（grad）
以文件名tidu_fuhao.m保存。
步骤二：
编写梯度函数。
打开Editor编辑器，输入以下语句：
functionf=f1_tidu（x0）
x=x0
（1）;
y=x0
（2）;
f=[8*x^3+8*x*y^2-8*x-4*x*y
8*x^2*y+16*y^3-14*y-2*x^2]';
步骤三：
编写最速下降法的方法函数打开Editor编辑器，输入以下语句：
function[x,n,data]=zuisu（x0）
ifnargin==1
tol=1e-6;
end
x1=x0-0.001*f1_tidu（x0）;
n=1;
%迭代过程
while（norm（x1-x0）>tol）&（n<2000）
x0=x1;
%0.001为步长因子
x1=x0-0.001*f1_tidu（x0）;
n=n+1;
%data用来存放中间数据
data（:
n）=x1;
end
x=x1;
以文件名zuisu.m保存。
步骤四：
编写主函数。
打开Editor编辑器，输入以下语句：
%最速下降法的主函数
x0=[0.80.7];
[x,n,data]=zuisu（x0）;
disp（'用最速下降法的计算结果为'）
x
disp（'迭代次数为'）
n
运行结果：
运行结果分析：
由两种算法的运行结果截屏可看出，同样的精度要求（10-5）下，Newton法收敛速度远远好于最速下降法，我们知道最速下降法是线性收敛的。
上机作业4
1、用幂法与反幂法求矩阵A的按模最大、最小特征值与对应的特征向量。
e=10-5
幂法原理
当矩阵A满足一定条件时，在工程中可用幂法计算其主特征值（按模最大）及其特征向量。
矩阵A需要满足的条件为：
（1）
（2）存在n个线性无关的特征向量，设为
不全为0，则有
可见，当
越小时，收敛越快；且当k充分大时，有
，对应的特征向量即是
。
设An有n个线性相关的特征向量v1，v2，…，vn，对应的特征值1，2，…，n，满足
|1|>|2|…|n|
（3.2.1）
规范化
在实际计算中，若|1|>1则|1ka1|，若|1|<1则|1ka1|0都将停机。
须采用“规范化”的方法
，k=0，1，2，…
定理3.2-1任给初始向量
有，
程序：
function[t,y]=lpower（A,x0,eps,N）
A=[4,1,1,1;1,3,-1,1;1,-1,2,0;1,1,0,2];
x0=[1;1;1;1];
eps=1e-5;N=100;k=1;z=0;
y=x0./max（abs（x0））;
x=A*y;
b=max（x）;
ifabs（z-b）t=max（x）;
return;
end
whileabs（z-b）>eps&&kk=k+1;
z=b;
y=x./max（abs（x））;
x=A*y;
b=max（x）;
end
[m,index]=max（abs（x））;
t=x（index）;
disp（'用幂法求矩阵的最大特征值及对应的特征向量'）
disp（'最大特征值为：
'）
t
disp（'最大特征值对应的特征向量为：
'）
y
disp（'迭代次数：
'）
k
end
反幂法理论
在工程计算中，可以利用反幂法计算矩阵按模最小特征值及其对应特征向量。
其基本理论如下，与幂法基本相同：
，可知，A和A-1的特征值互为倒数，求A按模最小特征值即求A-1的按模最大特征值,取倒数即为A的按模最小特征值所以算法基本相同，区别就是在计算
程序：
function[s,y]=invpower（A,x0,eps,n）
A=[4,1,1,1;1,3,-1,1;1,-1,2,0;1,1,0,2];
x0=[1;1;1;1];
eps=10^（-5）;
N=100;
k=1;
r=0;
y=x0./max（abs（x0））;
[L,U]=lu（A）;
z=L\y;
x=U\z;
u=max（x）;
s=1/u;
whileabs（u-r）>eps&&kk=k+1;
r=u;
y=x./max（abs（x））;
z=L\y;
x=U\z;
u=max（x）;
end
[m,index]=max（abs（x））;
s=1/x（index）;
disp（'用反幂法求矩阵的最小特征值及对应的特征向量'）
disp（'最小特征值为：
'）
s
disp（'最小特征值对应的特征向量为：
'）
y
disp（'迭代次数：
'）
k
end
运行结果：
运行结果分析：
由上图可知，幂法与反幂法求解矩阵的特征值与特征向量是很方便的。
2、用Householder变换求矩阵A的QR分解，并用QR方法做3次迭代。
原理：
对任意一个非奇异矩阵（可逆矩阵）A，可以把它分解成一个正交阵Q和一个上三角阵的乘积，称为对矩阵A的QR分解，即A=QR。
如果规定R的对角元取正实数，这种分解是唯一的。
若A是奇异的，则A有零特征值。
任取一个不等于A的特征值的实数μ，则A‐μI是非奇异的。
只要求出A‐μI的特征值和特征向量就容易求出矩阵A的特征值和特征向量，所以假设A是非奇异的，不是一般性。
设A=A1，对A1作QR分解，得A1=Q1R1，交换该乘积的次序，得A2=R1Q1=，由于Q1正交矩阵，A1到A2的变换为正交相似变换，于是A1和A2就有相同的特征值。
一般的令A1=A，对k=1，2，3，…..
这样，可得到一个迭代序列{Ak}，这就是QR方法的基本过程。
QR算法的核心思想是：
对一个具有n个不相等的特征值的矩阵A，进行如下变换：
A=Q1R1;A1=R1Q1
A1=Q2R2;A2=R2Q2
…………
Ak=Qk+1Rk+1;Ak+1=Rk+1Qk+1
…………
其中Qi是正交矩阵，Ri是可逆的上三角矩阵。
由于：
Ai+1=QTi+1（Qi+1Ri+1）Qi+1=QTi+1AiQi+1
所以Ai+1和Ai相似，具有相同的特
升级会员 