Jacobi迭代法GaussSeidel迭代法.docx
- 文档编号:8732275
- 上传时间:2023-02-01
- 格式:DOCX
- 页数:18
- 大小:125.69KB
Jacobi迭代法GaussSeidel迭代法.docx
《Jacobi迭代法GaussSeidel迭代法.docx》由会员分享,可在线阅读,更多相关《Jacobi迭代法GaussSeidel迭代法.docx(18页珍藏版)》请在冰豆网上搜索。
Jacobi迭代法GaussSeidel迭代法
Matlab线性方程组的迭代解法(Jacobi迭代法Gauss-Seidel迭代法)实验报告
2008年11月09日星期日12:
49
1.熟悉Jacobi迭代法,并编写Matlab程序matlab程序
按照算法(Jacobi迭代法)编写Matlab程序(Jacobi.m)
function[x,k,index]=Jacobi(A,b,ep,it_max)
%求解线性方程组的Jacobi迭代法,其中
%A---方程组的系数矩阵
%b---方程组的右端项
%ep---精度要求。
省缺为1e-5
%it_max---最大迭代次数,省缺为100
%x---方程组的解
%k---迭代次数
%index---index=1表示迭代收敛到指定要求;
%index=0表示迭代失败
ifnargin<4it_max=100;end
ifnargin<3ep=1e-5;end
n=length(A);k=0;
x=zeros(n,1);y=zeros(n,1);index=1;
while1
fori=1:
n
y(i)=b(i);
forj=1:
n
ifj~=i
y(i)=y(i)-A(i,j)*x(j);
end
end
ifabs(A(i,i))<1e-10|k==it_max
index=0;return;
end
y(i)=y(i)/A(i,i);
end
ifnorm(y-x,inf) break; end x=y;k=k+1; end 用Jacobi迭代法求方程组 的解。 输入: A=[430;33-1;0-14]; b=[24;30;-24]; [x,k,index]=Jacobi(A,b,1e-5,100) 输出: x= -2.9998 11.9987 -3.0001 k= 100 index= 0 2.熟悉Gauss-Seidel迭代法,并编写Matlab程序 function[v,sN,vChain]=gaussSeidel(A,b,x0,errorBound,maxSp) %Gauss-Seidel迭代法求解线性方程组 %A-系数矩阵b-右端向量x0-初始迭代点errorBound-近似精度maxSp-最大迭代次数 %v-近似解sN-迭代次数vChain-迭代过程的所有值 step=0; error=inf; s=size(A); D=zeros(s (1)); vChain=zeros(15,3);%最多能记录15次迭代次数 k=1; fx0=x0; fori=1: s (1) D(i,i)=A(i,i); end; L=-tril(A,-1); U=-triu(A,1); whileerror>=errorBound&step x0=inv(D)*(L+U)*x0+inv(D)*b; vChain(k,: )=x0'; k=k+1; error=norm(x0-fx0); fx0=x0; step=step+1; end v=x0; sN=step; 用Gauss-Seidel迭代法求解上题的线性方程组,取。 输入: A=[430;33-1;0-14]; b=[24;30;-24]; x0=[0;0;0]; [v,sN,vChain]=gaussSeidel(A,b,x0,0.00001,11) 输出: v= 0.6169 11.1962 -4.2056 sN= 11 vChain= 6.000010.0000-6.0000 -1.50002.0000-3.5000 4.500010.3333-5.5000 -1.75003.6667-3.4167 3.250010.6111-5.0833 -1.95835.0556-3.3472 2.208310.8426-4.7361 -2.13196.2130-3.2894 1.340311.0355-4.4468 -2.27667.1775-3.2411 0.616911.1962-4.2056 000 000 000 000 s 数值实验 数值实验要求: 数值实验报告内容: 要包含题目、算法公式、完整的程序、正确的数值结果和图形以及相应的误差分析。 在本课程网站上提交数值实验报告的电子文档。 一、为了逼近飞行中的野鸭的顶部轮廓曲线,已经沿着这条曲线选择了一组点。 见下表。 1.对这些数据构造三次自然样条插值函数,并画出得到的三次自然样条插值曲线; 2.对这些数据构造Lagrang插值多项式,并画出得到的Lagrang插值多项式曲线。 x0.91.31.92.12.63.03.94.44.75.06.0 f(x)1.31.51.852.12.62.72.42.152.052.12.25 x7.08.09.210.511.311.612.012.613.013.3 f(x)2.32.251.951.40.90.70.60.50.40.25 1.使用三次样条插值函数csape()求解。 解: 输入: >>x=[0.91.31.92.12.63.03.94.44.75.06.07.08.09.210.511.311.612.012.613.013.3]; >>y=[1.31.51.852.12.62.72.42.152.052.12.252.32.251.951.40.90.70.60.50.40.25]; >>pp=csape(x,y,'variational',[00]) 得到: pp= form: 'pp' breaks: [0.90001.30001.90002.10002.600033.90004.40004.700056789.200010.500011.300011.60001212.60001313.3000] coefs: [20x4double] pieces: 20 order: 4 dim: 1 再输入: >>pp.coefs 得到: ans= -0.247600.53961.3000 0.9469-0.29720.42081.5000 -2.95641.40731.08681.8500 -0.4466-0.36661.29492.1000 0.4451-1.03650.59342.6000 0.1742-0.5025-0.02222.7000 0.0781-0.0322-0.50342.4000 1.31420.0849-0.47712.1500 -1.58121.2676-0.07132.0500 0.0431-0.15550.26232.1000 -0.0047-0.02610.08082.2500 -0.0244-0.04010.01462.3000 0.0175-0.1135-0.13902.2500 -0.0127-0.0506-0.33581.9500 -0.0203-0.1002-0.53181.4000 1.2134-0.1490-0.73120.9000 -0.83930.9431-0.49290.7000 0.0364-0.0640-0.14130.6000 -0.44800.0014-0.17890.5000 0.5957-0.5361-0.39280.4000 因此,三次样条函数S(x)为 最后输入: >>xx=0: 0.01: 14;yy=interp1(x,y,xx,'csape'); >>plot(xx,yy);xlabel('x');ylabel('y'); 得到图形: 2.Lagrange插值算法: 1)输入N个节点数组; 2)定义初始值和; 3)利用公式 求N次插值基函数; 4)利用多项式求插值函数; 解: 输入: >>x=[0.9,1.3,1.9,2.1,2.6,3.0,3.9,4.4,4.7,5.0,6.0,7.0,8.0,9.2,10.5,11.3,11.6,12.0,12.6,13.0,13.3]; y=[1.3,1.5,1.85,2.1,2.6,2.7,2.4,2.15,2.05,2.1,2.25,2.3,2.25,1.95,1.4,0.9,0.7,0.6,0.5,0.4,0.25]; >>a=polyfit(x,y,length(x)-1); >>p=vpa(poly2sym(a),5) 得到数值结果: p=.30732e-10*x^20+.42770e-8*x^19-.27708e-6*x^18+.11098e-4*x^17-.30784e-3*x^16+.62785e-2*x^15-.97558e-1*x^14+1.1810*x^13-11.297*x^12+86.085*x^11-524.68*x^10+2558.0*x^9-9942.3*x^8+30586.*x^7-73622.*x^6+.13627e6*x^5-.18907e6*x^4+.18913e6*x^3-.12803e6*x^2+52170.*x-9593.4 再输入: >>y1=polyval(a,x); plot(x,y1);xlabel('x');ylabel('y'); 得到图形: 结果分析: 由以上两图形可以看出三次样条插值的图形较lagrange插值的图形要光滑的多。 因为样条函数插值不仅具有较好的收敛性和稳定性,而且其光滑性也较高,因此样条函数成为了重要的插值工具,其中应用较多的是三次样条插值。 二、对于问题 将h=0.025的Euler法,h=0.05的改进的Euler法和h=0.1的4阶经典的Runge-Kutta法在这些方法的公共节点0.1,0.2,0.3,0.4和0.5处进行比较。 精确解为: 。 1.Euler法 在x,y平面上微分方程的解在曲线上一点(x,y)的切线斜率等于函数的值。 该曲线的顶点设为p,再推进到p(),显然两个顶点p,p的坐标有以下关系 Matlab程序: function[x,y]=Euler(ydot,x0,y0,h,N) %ydot为一阶微分方程的函数 %x0,y0为初始条件 %h为区间步长 %N为区间个数 %x为Xn构成的向量,y为Yn构成的向量 x=zeros(1,N+1);y=zeros(1,N+1); x (1)=x0;y (1)=y0; forn=1: N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(ydot,x(n),y(n)); end 解: 取h=0.025,N=20,输入: >>ydot=inline('y-x^2+1','x','y'); [t,u]=Euler(ydot,0,0.5,0.025,20) 得到数值结果: t= Columns1through17 00.02500.05000.07500.10000.12500.15000.17500.20000.22500.25000.27500.30000.32500.35000.37500.4000 Columns18through21 0.42500.45000.47500.5000 u= Columns1through17 0.50000.53750.57590.61530.65550.69660.73870.78160.82530.87000.91550.96181.00891.05691.10571.15531.2056 Columns18through21 1.25681.30871.36131.4147 即采用Euler法得到: u(0.1)=0.6555,u(0.2)=0.8253,u(0.3)=1.0089,u(0.4)=1.2056,u(0.5)=1.4147 2.改进Euler方法 改进Euler公式. y=yn+h/2(f()+f(xn+1,+h*f()))即迭代公式为: Matlab程序: function[x,y]=Euler_ad(ydot,x0,y0,h,N) %改进Euler公式 %ydot为一阶微分方程的函数 %x0,y0为初始条件 %h为区间步长 %N为区间个数 %x为Xn构成的向量,y为Yn构成的向量 x=zeros(1,N+1);y=zeros(1,N+1); x (1)=x0;y (1)=y0; forn=1: N x(n+1)=x(n)+h; ybar=y(n)+h*feval(ydot,x(n),y(n)); y(n+1)=y(n)+h/2*(feval(ydot,x(n),y(n))+feval(ydot,x(n+1),ybar)); end 解: 取h=0.05,N=10,输入: >>ydot=inline('y-x^2+1','x','y'); [t,u]=Euler_ad(ydot,0,0.5,0.05,10) 得到数值结果: t= 00.05000.10000.15000.20000.25000.30000.35000.40000.45000.5000 u= 0.50000.57680.65730.74140.82910.92021.01471.11261.21361.31781.4250 即采用改进的Euler法得到: u(0.1)=0.6573,u(0.2)=0.8291,u(0.3)=1.0147,u(0.4)=1.2136,u(0.5)=1.4250 3.四阶Runge_Kutta法 求解公式为: Matlab程序: function[x,y]=Runge_Kutta4(ydot,x0,y0,h,N) %标准四阶Runge_Kutta方法 %ydot为一阶微分方程的函数 %x0,y0为初始条件 %h为区间步长 %N为区间个数 %x为Xn构成的向量,y为Yn构成的向量 x=zeros(1,N+1);y=zeros(1,N+1); x (1)=x0;y (1)=y0; forn=1: N x(n+1)=x(n)+h; k1=h*feval(ydot,x(n),y(n)); k2=h*feval(ydot,x(n)+1/2*h,y(n)+1/2*k1); k3=h*feval(ydot,x(n)+1/2*h,y(n)+1/2*k2); k4=h*feval(ydot,x(n)+h,y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4); end 解: 取h=0.1,N=5,输入: >>ydot=inline('y-x^2+1','x','y'); [t,u]=Runge_Kutta4(ydot,0,0.5,0.1,5) 得到数值结果: t= 00.10000.20000.30000.40000.5000 u= 0.50000.65740.82931.01511.21411.4256 结果比较: tEuler法改进Euler法四阶Runge_Kutta精确解 0.10.65550.65730.65740.6574 0.20.82530.82910.82930.8293 0.31.00891.01471.01511.0151 0.41.20561.21361.21411.2141 0.51.41471.42501.42561.4256 由以上结果可以看出改进的Euler法较Euler法计算精度有所提高,但还不是十分精确。 四阶Runge_Kutta法具有非常高的精度,事实上,在求解微分方程初值问题,四阶Runge_Kutta法是单步长中最优秀的方法,通常都用该方法求解实际问题。 三、 用Newton迭代法求方程的根时,分别取初始值,; 用Newton迭代法求方程时,分别取初始值,; 算法: (1)取初始点x0最大迭代次数N和精度要求ε,k=0. (2)如果f’(xk)=0,则停止计算;否则计算 Xk+1=xk-f(xk)/f’(xk) (3)若|xk+1-xk|<ε,则停止计算。 (4)若k=N,则停止计算;否则置k=k+1,转 (2)。 Matlab程序: function[x_star,index,it]=Newton(fun,x,ep,it_max) %求解非线性方程的Newton法 %fun(x)为需要求根的函数,第一个分量是函数值,第二个分量是导数值 %fun=inline('[x^3-x-1,3*x^2-1]')当f(x)=x^3-x-1; %x为初始点 %ep为精度,缺省值为1e-5 %it_max为最大迭代次数,缺省值100 %x_star为当迭代成功时输出方程的根,失败时输出最后的迭代值 %index为指标变量,index=1表示迭代成功index=0表示失败 %it为迭代次数 ifnargin<4it_max=100;end ifnargin<3ep=1e-5;end index=0;k=1; whilek<=it_max x1=x;f=feval(fun,x); ifabs(f (2)) x=x-f (1)/f (2); ifabs(x-x1) k=k+1; end x_star=x;it=k; 解: (1)由于f(x)=arctan(x),f’(x)=1/1+x2,取初始值x0=1.45,输入 >>fun=inline('[atan(x),1/(1+x^2)]'); >>[x_star,index,it]=Newton(fun,1.45) 得到数值结果: x_star=1.6281e+004 index=0 it=7 取初始值x0=0.5,输入 >>fun=inline('[atan(x),1/(1+x^2)]'); >>[x_star,index,it]=Newton(fun,0.5) 得到数值结果: x_star=0 index=1 it=4 输入x=-3: 0.05: 3;y=atan(x); plot(x,y);xlabel('x');ylabel('y'); 得到图形: (2)由于f(x)=x3-x-3=0,f’(x)=3x2-1,取初始值x0=0.0,输入 >>fun=inline('[x^3-x-3,3*x^2-1]'); >>[x_star,index,it]=Newton(fun,0.0) 得到数值结果: x_star=-0.0074 index=0 it=101 取初始值x0=2.0,输入 >>fun=inline('[x^3-x-3,3*x^2-1]'); >>[x_star,index,it]=Newton(fun,2.0) 得到数值结果: x_star=1.6717 index=1 it=4 输入x=0: 0.05: 3;y=x.^3-x-3; plot(x,y);xlabel('x');ylabel('y'); 得到图形: 结果分析: 从以上结果可以看出初值的选取与Newton法的收敛很有关系。 初值选的不好,Nexton法将不收敛,它的收敛性是在跟a附近讨论的。 选取初值时一定要十分小心。 四、用Jacobi迭代和SOR法分别求解线性方程组(教科书第86页算例2),并验证、输出SOR法的松弛因子w和对应的迭代次数。 1.Jacobi迭代法 Jacobi迭代法的算法为: Matlab程序: function[x,k,index]=Jacobi(A,b,ep,it_max) %求解线性方程组的Jacobi迭代法 %A为系数矩阵 %b为方程组右端项 %ep为精度要求,缺省值1e-5 %it_max为最大迭代次数,缺省值100 %x为方程组的解 %k为迭代次数 %index为指标变量index=1表示迭代收敛到指定要求index=0表示迭代失败。 ifnargin<4it_max=100;end ifnargin<3ep=1e-5;end n=length(A);k=0; x=zeros(n,1);y=zeros(n,1);index=1; while1 fori=1: n y(i)=b(i); forj=1: n ifj~=i y(i)=y(i)-A(i,j)*x(j); end end ifabs(A(i,i))<1e-10|k==it_max index=0;return; end y(i)=y(i)/A(i,i); end ifnorm(y-x,inf) break; end x=y;k=k+1; end function[A,b]=shuru %自己定义输入矩阵A和向量b的函数 n=15; fori=1: n ifi==1|i==15z(i)=3; elsez(i)=2; end forj=1: n ifj
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Jacobi 迭代法 GaussSeidel