计算方法MATLAB实验报告Runge现象与方程求根合集

实验Runge现象的演示【实验目的】通过上机绘制Runge函数图像，理解高次插值的病态性质。【实验内容】在区间[-1,1]上取n=11等距节点对龙格(Runge)函数f(x)=1/(1+25*x^2)做多项式插值，并画出插值函数。【实验程序】functiontry_Rungef=inline('1./(1+25*x.^2)');n=11;X=linspace(-1,1,n);Y=f(X);x=-1:0.01:1;y=mylagpoly(X,Y,x);plot(x,f(x),'r',X,Y,'o',x,y,'b')title('Runge现象')legend('y=1/(1+25*x^2)','插值点','等分的10次插值多项式',0)function[C,L]=interp_lag(X,Y)n=length(X);L=zeros(n,n);forj=1:nC1=[X(1:j-1),X(j+1:n)];C2=poly(C1);A=polyval(C2,X(j));L(j,:)=C2/A;endC=Y*L;functiony=mylagpoly(X,Y,x)n=length(X);y=zeros(size(x));fori=1:nt=1;forj=1:nifj~=it=t.*(x-X(j))/(X(i)-X(j));endendy=y+Y(i)*t;end【运行结果如下】：【结果分析】：在实际运用中，为了取得较好的插值逼近效果，应尽量保证以下几点：不采用次数过高的插值多项式；适当选取插值节点；避免函数值突变，若不得已对存在不连续点的函数进行插值逼近，可以尝试分段插值，并将不连续点都处理到子区间的端点上，从而原函数在各子区间内分段连续，以便提高插值逼近效果。

专业机械序号姓名日期2017-12-20实验二非线性方程求根【实验目的】初步认识非线性问题的迭代法与线性问题迭代法的差别，探讨迭代法及初始值与迭代收敛性的关系。【问题提出】迭代法是求解非线性方程的基本思想方法，与线性方程的情况一样，其构造方法可以有多种多样，但关键是怎样才能使迭代收敛且有较快的收敛速度。【实验内容】考虑一个简单的代数方程针对上述方程，可以构造多种迭代法，如在实轴上取初始值x0，请分别用迭代（7.1）-（7.3）作实验，记录各算法的迭代过程。实验要求：（1）取定某个初始值，分别计算（7.1）-（7.3）迭代结果，它们的收敛性如何？重复选取不同的初始值，反复实验。请自选设计一种比较形象的记录方式（如利用MATLAB的图形功能），分析三种迭代法的收敛性与初值选取的关系。（2）对三个迭代法中的某个，取不同的初始值进行迭代，结果如何？试分析迭代法对不同的初值是否有差异？2.教材第二章课后实验课题实验要求：取定某个初始值，比较牛顿法和其他迭代法的数值结果，它们的收敛性如何？重复选取不同的初始值，反复实验。请自选设计一种比较形象的记录方式（如利用MATLAB的图形功能），分析三种迭代法的收敛性与初值选取的关系。【结果及分析】:实验1程序如下：functiontry_stableglobalnN=20;%计算N个值a=1.6;%初值%%--------------------------------------------%%[方案I]用递推公式I0=a;%初值I=zeros(N,1);%创建Nx1矩阵(即列向量),元素全为零I(1)=I0^2-1;forn=1:NI(n+1)=I(n)^2-1;err=abs(I(n+1)-I(n));%相邻两次迭代的误差if(err<0.000001)%'|'是'或'breakendend%%--------------------------------------------%%[方案II]用递推公式II0=a;%初值II=zeros(N,1);%创建Nx1矩阵(即列向量),元素全为零II(1)=1+1/II0;forn=1:NII(n+1)=1+1/II(n);err=abs(II(n+1)-II(n));%相邻两次迭代的误差if(err<0.000001)%'|'是'或'breakendend%%--------------------------------------------%%调用matlab高精度数值积分命令quadl计算以便比较%III0=a;%初值III=zeros(N,1);%创建Nx1矩阵(即列向量),元素全为零III(1)=(1+III0)^0.5;forn=1:NIII(n+1)=(1+III(n))^0.5;err=abs(III(n+1)-III(n));%相邻两次迭代的误差if(err<0.000001)%'|'是'或'breakendendf1=inline('x^2-x-1');X2=fzero(f1,[12])%%--------------------------------------------%%显示计算结果%clcfprintf('\n方案I结果方案II结果方案III结果精确值')forn=1:N,fprintf('\nI(%2.0f)%17.7f%17.7f%17.7f%17.7f',n,I(n),II(n),III(n),X2)end运行结果如下：方案I结果方案II结果方案III结果精确值I(1)1.56000001.62500001.61245151.6180340I(2)1.43360001.61538461.61630801.6180340I(3)1.05520901.61904761.61750051.6180340I(4)0.11346591.61764711.6178691.6180340I(5)-0.98712551.61818181.61798301.6180340I(6)-0.0255831.61797751.61801821.6180340I(7)-0.99934551.61805561.61802911.6180340I(8)-0.00130861.61802581.61803251.6180340I(9)-0.99999831.61803711.61803351.6180340I(10)-0.00000341.61803281.61803381.6180340I(11)-1.00000001.61803440.00000001.6180340I(12)-0.00000001.61803380.00000001.6180340I(13)-1.00000000.00000000.00000001.6180340I(14)0.00000000.00000000.00000001.6180340I(15)-1.00000000.00000000.00000001.6180340I(16)0.00000000.00000000.00000001.6180340I(17)-1.00000000.00000000.00000001.6180340I(18)0.00000000.00000000.00000001.6180340I(19)-1.00000000.00000000.00000001.6180340I(20)0.0000000.00000000.00000001.6180340>>结果分析：方案一，开始误差较小，误差逐渐增大后开始振荡，说明方案不可行。方案二，开始误差较小，但直到第11步才接近真值，说明其快速性不够好。方案三，开始误差较小，较于方案二更快速。实验2程序如下：functiontry_stableglobalnN=20;%计算N个值a=3;%%--------------------------------------------%%[方案I]用牛顿法%I0=a;%初值I=zeros(N,1);%创建Nx1矩阵(即列向量),元素全为零I(1)=I0-(I0^3-sin(I0)-12*I0+1)/(3*I0^2-cos(I0)-12);forn=1:NI(n+1)=I(n)-(I(n)^3-sin(I(n))-12*I(n)+1)/(3*I(n)^2-cos(I(n))-12);err=abs(I(n+1)-I(n));%相邻两次迭代的误差if(err<0.000001)breakendend%%--------------------------------------------%%[方案II]用递推公式%%迭代公式一x=φ_1(x)=√(3&12x+sinx-1)II0=a;%初值II=zeros(N,1);%创建Nx1矩阵(即列向量),元素全为零II(1)=(12*II0+sin(II0)-1)^(1/3);forn=1:NII(n+1)=(12*II(n)+sin(II(n))-1)^(1/3);err=abs(II(n+1)-II(n));%相邻两次迭代的误差if(err<0.000001)breakendend%%--------------------------------------------%%迭代公式二x=φ_2(x)=1/12(x^3-sinx+1)III0=a;%初值III=zeros(N,1);%创建Nx1矩阵(即列向量),元素全为零III(1)=(1/12)*(III0^3-sin(III0)+1);forn=1:NIII(n+1)=(1/12)*(III(n)^3-sin(III(n))+1);err=abs(III(n+1)-III(n));%相邻两次迭代的误差if(err<0.000001)breakendend%%[方案III]用埃特金迭代加速法IIII0=a;%初值IIII=zeros(N,1);%创建Nx1矩阵(即列向量),元素全为零IIII(1)=(12*IIII0+sin(IIII0)-1)^(1/3);IIII(2)=(12*IIII(1)+sin(IIII(1))-1)^(1/3)forn=1:NIIII(n+1)=(12*IIII(n)+sin(IIII(n))-1)^(1/3);IIII(n+2)=(12*IIII(n+1)+sin(IIII(n+1))-1)^(1/3);IIII(n+1)=IIII(n+2)-((IIII(n+2)-IIII(n+1))^2)/(IIII(n+2)-2*IIII(n+1)+IIII(n));err=abs(IIII(n+1)-IIII(n));%相邻两次迭代的误差if(err<0.000001)breakendend%%--------------------------------------------f1=inline('x^2-x-1');X2=fzero(f1,[12])%%--------------------------------------------%%显示计算结果%clcfprintf('\n牛顿法结果普通法(1)结果普通法(2)结果埃特金加速法结果')forn=1:N,fprintf('\nI(%2.0f)%17.7f%17.7f%17.7f%17.7f',n,I(n),II(n),III(n),IIII(n))end运行结果如下：牛顿法结果普通法(1)结果普通法(2)结果埃特金加速法结果I(1)3.50913853.27545672.32157333.2754567I(2)3.41401403.36704071.06511933.4103566I(3)3.41013153.39644800.11112573.4101251I(4)3.41012513.40579380.07420633.4101251I(5)3.41012513.40875450.07718923.4101251I(6)0.00000003.40969150.07694560.0000000I(7)0.00000003.40998790.07696550.0000000I(8)0.00000003.41008170.07696390.0000000I(9)0.0000003.41011130.07696400.0000000I(10)0.00000003.4101200.00000000.0000000I(11)0.00000003.41012370.00000000.0000000I(12)0.00000003.41012460.00000000.0000000I(13)0.00000000.00000000.00000000.0000000I(14)0.00000000.00000000.000000.0000000I(15)0.00000000.00000000.00000000.0000000I(16)0.00000000.00000000.00000000.0000000I(17)0.00000000.00000000.00000000.0000000I(18)0.00000000.00000000.00000000.0000000I(19)0.00000000.00000000.00000000.0000000I(20)0.00000000.00000000.00000000.0000000>>【结果分析】牛顿法，开始误差较小收敛性好，迭代速度快，跟艾特金加速法几乎相同。但开始误差不稳定。普通方法1，快速性较差，结果存在一定误差。普通方法2，误差较大，还被进一步放大，不可取。

专业机械序号姓名日期2017-12-20实验三病态的线性方程组的求解【实验目的】了解病态方程组的特点，编程用迭代格式求解Hilbert矩阵掌握高斯、雅可比、高斯—赛德尔三种迭代格式；掌握用如何用MATLAB程序表示这三种迭代格式，并用其求解问题。【实验内容】问题提出：理论的分析表明，求解病态的线性方程组是困难的。实际情况是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b的求解，其中系数矩阵H为Hilbert矩阵，这是一个著名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端b的办法给出确定的问题。实验要求：（1）选择问题的维数为6，分别用Gauss消去法、Jacobi迭代法、GS迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？（2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？【程序】,求解病态方程组是相当困难%一般的方法都会失效,例如列主元Gauss消去法,要采取特殊的方法和技术%例如pcg(预处理共轭梯度法)等方法%Hilbert矩阵是严重的病态的矩阵(见P56)n=6;%矩阵阶数（迭代次数）A=hilb(n);%n阶Hilbert矩阵x=ones(n,1);%指定精确解全是1b=A*x;c=cond(A,inf);%求条件数x1=A\b;%用Gauss列主元消去法求解%%%%%%%%%雅可比V=diag(A);D=diag(V);L=-tril(A,-1);U=-triu(A,1);x2=zeros(n,1);%初值fork=0:100x2=inv(D)*(L+U)*x2+inv(D)*b;end%%%%%%%%%%%%高斯-塞德尔x3=zeros(n,1);fori=0:100x3=inv(D-L)*U*x3+inv(D-L)*bend%显示结果clcfprintf('A的条件数=%f\n',c)fprintf('\nGauss法求解结果雅可比法求解结果')fori=1:nfprintf('\n%6.2f%6.2f',x1(i),x2(i))endfprintf('\n高斯-塞德