Jacobi迭代法 Gauss-Seidel迭代法

Matlab线性方程组的迭代解法（Jacobi迭代法Gauss-Seidel迭代法）实验报告2008年11月09日星期日12:491.熟悉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表示迭代失败ifnargin4it_max=100;endifnargin3ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while1fori=1:ny(i)=b(i);forj=1:nifj~=iy(i)=y(i)-A(i,j)*x(j);endendifabs(A(i,i))1e-10|k==it_maxindex=0;return;endy(i)=y(i)/A(i,i);endifnorm(y-x,inf)epbreak;endx=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.999811.9987-3.0001k=100index=02.熟悉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&stepmaxSpx0=inv(D)*(L+U)*x0+inv(D)*b;vChain(k,:)=x0';k=k+1;error=norm(x0-fx0);fx0=x0;step=step+1;endv=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.616911.1962-4.2056sN=11vChain=6.000010.0000-6.0000-1.50002.0000-3.50004.500010.3333-5.5000-1.75003.6667-3.41673.250010.6111-5.0833-1.95835.0556-3.34722.208310.8426-4.7361-2.13196.2130-3.28941.340311.0355-4.4468-2.27667.1775-3.24110.616911.1962-4.2056000000000000s数值实验数值实验要求：数值实验报告内容：要包含题目、算法公式、完整的程序、正确的数值结果和图形以及相应的误差分析。在本课程网站上提交数值实验报告的电子文档。一、为了逼近飞行中的野鸭的顶部轮廓曲线，已经沿着这条曲线选择了一组点。见下表。1．对这些数据构造三次自然样条插值函数，并画出得到的三次自然样条插值曲线；2．对这些数据构造Lagrang插值多项式，并画出得到的Lagrang插值多项式曲线。x0.91.31.92.12.63.03.94.44.75.06.0f(x)1.31.51.852.12.62.72.42.152.052.12.25x7.08.09.210.511.311.612.012.613.013.3f(x)2.32.251.951.40.90.70.60.50.40.251．使用三次样条插值函数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:20order:4dim:1再输入：pp.coefs得到：ans=-0.247600.53961.30000.9469-0.29720.42081.5000-2.95641.40731.08681.8500-0.4466-0.36661.29492.10000.4451-1.03650.59342.60000.1742-0.5025-0.02222.70000.0781-0.0322-0.50342.40001.31420.0849-0.47712.1500-1.58121.2676-0.07132.05000.0431-0.15550.26232.1000-0.0047-0.02610.08082.2500-0.0244-0.04010.01462.30000.0175-0.1135-0.13902.2500-0.0127-0.0506-0.33581.9500-0.0203-0.1002-0.53181.40001.2134-0.1490-0.73120.9000-0.83930.9431-0.49290.70000.0364-0.0640-0.14130.6000-0.44800.0014-0.17890.50000.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:Nx(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=Columns1through1700.02500.05000.07500.10000.12500.15000.17500.20000.22500.25000.27500.30000.32500.35000.37500.4000Columns18through210.42500.45000.47500.5000u=Columns1through170.50000.53750.57590.61530.65550.69660.73870.78160.82530.87000.91550.96181.00891.05691.10571.15531.2056Columns18through211.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.41472.改进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:Nx(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)得到数值结果