Jacobi迭代法.docx
《Jacobi迭代法.docx》由会员分享，可在线阅读，更多相关《Jacobi迭代法.docx（14页珍藏版）》请在冰豆网上搜索。
Jacobi迭代法
第一题
算法解释
Jacobi迭代法
方程组Ax=b，其中A∈Rnxn，b∈Rn，且A为非奇异，则A可以写成A=D-L-U。
其中，D=diag[a11,a22,…,ann]，而—L，—U分别为A的上三角和下三角部分（不包括对角线元素）
则x=D—1（L+U）x+D-1b，由此可以构造迭代法:
x（k+1）=Bx（k）+f
其中：
B=D-1（L+U）x=I-D-1A，f=D-1b.
M文件
function［x,n］=jacobi（A，b，x0,eps,varargin）
%采用Jacobi迭代法求线性方程组Ax=b的解
％线性方程组的系数矩阵：
A
%线性方程组中的常数向量：
b
％迭代初始向量：
x0
%解的精度控制：
eps
%迭代步数控制:
varargin
%线性方程组的解：
x
%求出所需精度的解实际的迭代步数:
n
ifnargin==3
eps=1.0e—6;
M=200;
elseifnargin〈3
error
return
elseifnargin==5
M=varargin｛1};
end
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；％迭代次数
％迭代过程
whilenorm（x—x0）>=eps
x0=x;
x=B＊x0+f;
n=n+1;
if（n>=M）
disp（’Warning:
迭代次数太多，可能不收敛！
’）；
return；
end
end
x
n
例题
Jacobi迭代法求线性方程组实例。
用Jacobi迭代法求解以下线性方程组
10x1—x2=9
—x1+10x2+—2x3=7
-x2+10x3=6
在matlab命令窗口输入如下程序：
〉>a=［10-10；—110-2;0—210］；
〉>b=[9;7；6]；
〉>jacobi（a，b，［0;0；0]）
y=
0.9958
0.9579
0。
7916
输出的迭代次数为：
n=
11
ans=
0.9958
0.9579
0.7916
n=length（x）
第二题
方法一：
运用戴维宁定理求出外电路的等效电压Uoc=（Z2/（Z1+Z2）-Z4/（Z3+Z4））.＊Us，将电压源短路，z1,z2并联，z3，z4并联，Req=Z3。
＊Z4./（Z3+Z4）+Z1。
*Z2./（Z1+Z2）。
所求电流源两端电压为Ut=Is。
＊Req+Uoc
w=［eps，1，2]；Us=[10,10,0］;Is=[5，5,0]；
Z1=1。
/（0.5*w＊j）;Z4=1*w*j;
Z2=［2,2，2］；Z3=[2,2,2];
Uoc=（Z2/（Z1+Z2）-Z4/（Z3+Z4））.*Us；
Req=Z3.＊Z4./（Z3+Z4）+Z1.*Z2./（Z1+Z2）;
Ut=Is。
＊Req+Uoc
disp（‘wUmphi‘）
disp（［w'，abs（u’）,angle（u’）＊180/pi]）
等效电路图
方法二:
运用电路叠加定理，分别求出电流源与电压源单独作用在所求端点两端的电压，当只有电压源作用时，电流源断路，得到Uoc=（z2./（z1+z2）—z4。
/（z3+z4））。
＊us;当只有电流源作用时，电压源短路,z1,z2并联,z3，z4并联,外加电阻Ｒ＝（z3。
*z4./（z3+z4）+z1.＊z2。
/（z1+z2），电流源作用得到Uc=Is.*Ｒ。
所求端点电压为两者叠加值Ut＝Ｕｏｃ＋Ｕｃ
〉>w=［eps，1，2］； ％输入频率值
〉〉us=[10，10，0]; %输入不同频率下的电压分量
>〉Is=［5，0,5］；
〉〉z1=1./（0。
5＊w*j）;
>>z4=1＊w*j;
>〉z2=[2,2,2］;
〉〉z3=［2，2,2];
>〉Uoc=（z2。
/（z1+z2）-z4。
/（z3+z4））。
＊us； %电压源单独作用的端点电压
>>Uc=Is.＊（z3.*z4。
/（z3+z4）+z1.＊z2./（z1+z2））；％电流源单独作用的端点电压
〉〉Ut=Uoc+Uc；
〉〉disp（[w’，abs（u'），angle（u’）＊180/pi］）
0。
000010.00000
1.00003。
1623—18。
4349
2.00007.0711-8.1301
开始
第三题:
已知f（x）=3sin（pi*x/6）.^2,x=1。
5，3.5
k
Xk
f（Xk）
0
1
2
3
4
0。
0
1.0
2.0
3。
0
4.0
0.000
0。
75
2.25
3.0
2.25
求：
（1）计算函数差商表。
（2）写出牛顿多项式Nk（x）,k=1，2，3,4.
（3）在给定值x处求牛顿多项式Nk（x）,k=1,2,3,4的值。
（4）比较（3）中的结果与实际函数值。
M文件
function[P,A]=chashang（X，Y）
n=length（X）；
A=zeros（n,n）；
A（:
，1）=Y’;
forj=2:
n
fori=j:
n
A（i，j）=（A（i,j-1）—A（i—1，j-1））/（X（i）—X（i—j+1））;
end
end
P=A（1:
n,1：
n）
在命令窗口中输入
>〉X=[01234］;
〉〉Y=[00.752。
253.02.25];
>〉chashang（X，Y）
P=
00000
0。
75000。
7500000
2。
25001。
50000。
375000
3.00000。
7500-0.3750-0.25000
2。
2500-0.7500—0.7500—0。
12500。
0313
ans=
00000
0.75000。
7500000
2。
25001.50000.375000
3.00000.7500—0。
3750-0。
25000
2.2500—0.7500-0.7500—0。
12500.0313
差商程序流程图
M文件
function[C,A]=newtonpoly（X，Y）
n=length（X）；
A=zeros（n,n）；
A（：
，1）=Y'；
forj=2:
n
fori=j：
n
A（i,j）=（A（i，j-1）-A（i—1,j—1））/（X（i）-X（i—j+1））；
end
end
C=A（n,n）;
fork=（n—1）:
—1：
1
C=conv（C,poly（X（k）））;
d=length（C）；
C（d）=C（d）+A（k，k）；
end
在命令窗口中输入
〉>X=［01］；
〉〉Y=[00。
75]；
>>n1=newtonpoly（X，Y）
n1=
0。
75000
〉〉N1=poly2sym（n1）
N1=
3/4*x
>〉X=[012］；
>〉Y=［00.752。
25］；
>>n2=newtonpoly（X，Y）
n2=
0。
37500。
37500
>〉N2=poly2sym（n2）
N2=
3/8＊x^2+3/8*x
>>X=［0123]；
>〉Y=［00.752.253.0］；
>>n3=newtonpoly（X，Y）
n3=
—0。
25001.1250—0。
12500
〉〉N3=poly2sym（n3）
N3=
—1/4*x^3+9/8＊x^2-1/8*x
>>X=[01234];
>〉Y=[00。
752。
253。
02.25];
〉〉n4=newtonpoly（X，Y）
n4=
0.0313-0。
43751.4688-0。
31250
〉〉N4=poly2sym（n4）
N4=
1/32*x^4-7/16＊x^3+47/32＊x^2—5/16＊x
>〉f1=polyval（n1,1.5）
f1=
1.1250
〉〉f2=polyval（n2，1。
5）
f2=
1.4063
>〉f3=polyval（n3,1.5）
f3=
1。
5000
>〉f4=polyval（n4,1。
5）
f4=
1.5176
〉>p1=polyval（n1,3.5）
p1=
2。
6250
>〉p2=polyval（n2,3.5）
p2=
5.9063
〉〉p3=polyval（n3，3。
5）
p3=
2.6250
〉>p4=polyval（n4，3.5）
p4=
2.8301
>〉x=1。
5；
〉>F=3*（sin（pi＊x/6））^2
F=
1.5000
〉>x=3.5；
〉>P=3*（sin（pi＊x/6））^2
P=
2.7990
画图比较
输入:
x=0：
0。
05:
4；
y=3.*sin（pi.*x。
/6）。
＊sin（pi.*x./6）;
>〉plot（x,y，’r’）
>〉holdon
>>x=1.5；
>〉f1=1。
125；
>>f2=1.4063
>>f3=1。
500;
〉〉f4=1.5176
〉〉plot（x，f1，’o’，x,f2，’o',x，f3，'o',x,f4，’o'）
〉〉x=3。
5；
>>p1=2。
6250；
〉〉p2=5。
9063；
〉〉p3=2。
6250;
〉>p4=2.8301；
>〉plot（x，p1，’o’,x,p2,'o’，x，p3,'o',x，p4,'o'）
得到如下图形
升级会员 