数值分析代码详解：多项式插值、曲线拟合、数值积分、线性方程组求解

数值分析代码详解：多项式插值、曲线拟合、数值积分、线性方程组求解/n/n本文提供了一系列数值分析代码，涵盖多项式插值、曲线拟合、数值积分以及线性方程组的直接解法和迭代解法。详细解释每个代码段的功能和实现原理，并附带图形展示结果。/n/n### 多项式插值/n/n第1题：/n/nmatlab/nclear;/nclc;/nclf;/nx1=[0.2 0.4 0.6 0.8 1.0];/ny1=[0.98 0.92 0.81 0.64 0.38];/nn=length(y1);/nc=y1(:);/nfor j=2:n %求差商/n     for i=n:-1:j/n        c(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));/n    end/nend/nsyms x df d;/ndf(1)=1;d(1)=y1(1);/nfor i=2:n %求牛顿差值多项式/n    df(i)=df(i-1)*(x-x1(i-1));/n    d(i)=c(i)*df(i);/nend/nP4=vpa(sum(d),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数/npp=csape(x1,y1, 'variational');%调用三次样条函数/nq=pp.coefs;/nfor i=1:4/n    S=q(i,:)*[(x-x1(i))^3;(x-x1(i))^2;(x-x1(i))^1;(x-x1(i))^0];/n    S=vpa(collect(S),5)/nend/n/nfigure/nezplot(P4, [0.2,1.08]);/nhold on;/nx2=0.2:0.08:1.08;/ny2=fnval(pp, x2);/nind=[1,2,11,10];/ny3=fnval(pp,x2(ind));/nplot(x2,y2,'r',x2(ind),y3,'go')/ntitle('Newton interpolation with deg=4 and cubic splines');/nhold off;/n/n/n该代码使用牛顿插值法和三次样条函数对给定数据点进行插值。/n/n* 前半部分代码计算牛顿插值多项式系数，使用嵌套循环计算差商，最终得到4次牛顿插值多项式P4。/n* 后半部分代码使用csape函数拟合三次样条曲线，并绘制P4和三次样条曲线的图形。/n/n第2题：/n/nmatlab/nfunction y =lagrange (x0,y0,x)/n  m=length(x);/n  n=length(x0);/n  for i=1:m/n      z=x(i);/n      s=0;/n      for k=1:n/n          p=1.0;/n          for j=1:n/n              if j~=k/n                  p=p*(z-x0(j))/(x0(k)-x0(j));/n              end/n          end/n          s=p*y0(k)+s;/n      end/n      y(i)=s;/n  end/n%%%%%%%%%%%%%%%%%%%%/nclear;/nclc;/nclf;/nx0=[-1:0.02:1];/ny0=1./(1+25*x0.^2);/nplot(x0,y0,'b') %绘制原曲线/n/nhold on/nx1=linspace(-1,1,11);/ny1=1./(1+25*x1.^2);/ny0=lagrange(x1,y1,x0);/nplot(x0,y0,'--r') %插值曲线 /n/nx1=linspace(-1,1,21);/ny1=1./(1+25*x1.^2);/ny0=lagrange(x1,y1,x0);/nplot(x0,y0,'--g') %插值曲线 /n/n/n该代码使用拉格朗日插值法对给定数据点进行插值。/n/n* lagrange函数实现了拉格朗日插值算法。/n* 后半部分代码分别使用11个和21个数据点进行插值，并绘制出插值曲线和原曲线的图形。/n/n### 曲线拟合/n/n参考答案：/n/n第1题：/n/nmatlab/nclear;/nclc;/nclf;/nx = -1:0.01:1;/nxf = -1:0.2:1;/ny = 1./(1+25*x.^2);/nyf = 1./(1+25*xf.^2);/ny1 = polyfit(xf,yf,3);/nY =polyval(y1,x);/ny2 =polyval(y1,xf);/nplot(x,y,x,Y,'m', xf,yf,'or',xf,y2,'*b');% 输出原函数曲线以及拟合多项式曲线./nf = poly2str(y1,'x')% 将三次多项式拟合后得到的多项式的系数向量表示成对应的多项式的习惯表达式./ntitle('Curve Fitting with d=3','FontName','New Times Roman','FontSize',12);/nxlabel('x-axis','FontName','New Times Roman','FontSize',12);/nylabel('y-axis','FontName','New Times Roman','FontSize',12);/nlegend('Given curve','Fitting curve')/n/n/n该代码使用polyfit函数对给定数据点进行三次多项式拟合，并绘制出拟合曲线和原曲线的图形。/n/n* polyfit函数计算三次多项式的系数，polyval函数计算拟合曲线上的点。/n* 最后绘制出原曲线、拟合曲线和拟合点。/n/n第2题：/n/nmatlab/nclear;/nclf;/nclc;/nx=[0 0.1 0.2 0.3 0.5 0.8 1];/ny=[1 0.41 0.5 0.61 0.91 2.02 2.46];/np1=polyfit(x,y,3);/np2=polyfit(x,y,4);/nX = 0:0.01:1;/ny1=polyval(p1,X);/ny2=polyval(p2,X);/npoly2str(p1,'x')/npoly2str(p2,'x')/nplot(x,y,'ko',X,y1,'b-',X,y2,'r-')/nhold on/np3=polyfit(x,y,2);%观察图形与抛物线接近，故采用2次曲线拟合/ny3=polyval(p3,X);/npoly2str(p3,'x')/nplot(X,y3,'m-')/ntitle('Curve Fitting with different degree','FontName','New Times Roman','FontSize',12);/nxlabel('x-axis','FontName','New Times Roman','FontSize',12);/nylabel('y-axis','FontName','New Times Roman','FontSize',12);/nlegend('Fitted pts', 'd=3','d=4','d=2')/nhold off/n/n/n该代码使用polyfit函数对给定数据点分别进行三次、四次和二次多项式拟合，并绘制出拟合曲线的图形。/n/n* 代码首先使用polyfit分别计算三次、四次和二次多项式的系数，然后使用polyval计算拟合曲线上的点。/n* 最后绘制出拟合点和不同次数的拟合曲线。/n/n### 数值积分/n/nmatlab/n% 1. 定义函数/nfunction y=fun41(x)/ny=sqrt(x).*log(x);/nend/n/n% 2. 用复合梯形公式和复合辛普森公式计算数值积分/nclear;/nclc;/nh=0.001;     %h为步长,可分别令h=1,0.1,0.01,0.001/nn=1/h; /nt=0; /nfor i=1:n-1/n    t=t+fun41(i*h);/nend/nT=h/2*(0+2*t+fun41(1));/nT=vpa(T,10) /n% 以上为复合梯形公式/n/n% 以下为复合辛普森公式/ns1=0;/ns2=0;/nfor i=0:n-1/n    s1=s1+fun41(h/2+i*h);/nend/nfor i=1:n-1/n    s2=s2+fun41(i*h);/nend/nS=h/6*(0+4*s1+2*s2+fun41(1));/nS=vpa(S,10) /n/n% 3. 使用龙贝格算法计算数值积分/nclear;/nclc;/nm=16;/nh=1;/nT(1)=(0+fun41(1))*h/2/nfor i=2:m/n    h=h/2;/n    n=1/h;/n    t=0;/n    for j=1:2:n-1/n        t=t+fun41(j*h);/n    end/n    T(i)=T(i-1)/2+h*t;%梯形公式/nend/nfor i=1:m-1/n    for j=m:i+1/n        T(j)=4^i/(4^i-1)*T(j)-1/(4^i-1)*T(j-1);/n        %通过不断的迭代求得T(j)，即T表的对角线上的元素。/n    end/nend/nvpa(T(m),10)/n/n/n该代码使用复合梯形公式、复合辛普森公式和龙贝格算法分别计算给定函数fun41的数值积分。/n/n* 代码首先定义函数fun41。/n* 然后使用复合梯形公式和复合辛普森公式分别计算数值积分。/n* 最后使用龙贝格算法对数值积分进行迭代计算，得到更加精确的结果。/n/n### 直接求解线性方程组/n/n参考答案：/n/n1. 第一问程序：/n/nmatlab/nclear;/nclc;/nA=[10  -7        0 1/n   -3   2.099999 6 2/n    5  -1        5 -1/n    2   1        0 2];/nb=[8;5.900001;5;1];/n[m,n]=size(A); /nL=eye(n); /nU=zeros(n); /nflag='ok'; /nfor i=1:n /n    U(1,i)=A(1,i); /nend /nfor r=2:n /n   L(r,1)=A(r,1)/U(1,1); /nend /nfor i=2:n /n     for j=i:n /n         z=L(i,1:i-1)*U(1:i-1,j);/n         U(i,j)=A(i,j)-z; /n     end /n     if abs(U(i,i))<eps /n         flag='failure' /n         return; /n     end /n     for k=i+1:n /n         m=L(k,1:i-1)*U(1:i-1,i);/n         L(k,i)=(A(k,i)-m)/U(i,i); /n     end /nend /nL/nU/ny=L/b;x=U/y/ndetA=det(L)*det(U)/n/n/n/n该代码使用LU分解法求解线性方程组。/n/n* 代码首先对系数矩阵A进行LU分解，得到下三角矩阵L和上三角矩阵U。/n* 然后分别解Ly=b和Ux=y，得到线性方程组的解x。/n* 最后计算A的行列式值。/n/n第二问程序：/n/nmatlab/nclear;/nclc;/nA=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];/nb=[8;5.900001;5;1];/n[n,n] = size(A); /nx = zeros(n,1); /nP = [1:n];/nAug = [A,b]; %增广矩阵/n /nfor k = 1:n-1 /n    [piv,r] = max(abs(Aug(k:n,k))); %找列主元所在子矩阵的行r/n    r = r + k - 1; % 列主元所在大矩阵的行/n    if r>k /n        Aug([k,r],:)=Aug([r,k],:); /n        P([k,r])=P([r,k]);/n    end/n    if Aug(k,k)==0 /n        error('对角元出现0');/n    end/n    % 把增广矩阵消元成为上三角/n    for p = k+1:n /n        Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); /n    end/nend/n /n% 解上三角方程组/nA = Aug(:,1:n); /nb = Aug(:,n+1); /nx(n) = b(n)/A(n,n); /nfor k = n-1:-1:1 /n    x(k) = (b(k)-A(k,n:-1:k+1)*x(n:-1:k+1))/A(k,k);/nend/nP/nx/ndetA=det(A)/n/n/n该代码使用列主元高斯消元法求解线性方程组。/n/n* 代码首先将系数矩阵A和常数向量b合并为增广矩阵Aug。/n* 然后进行列主元高斯消元，将Aug消元为上三角矩阵。/n* 最后通过回代求解上三角矩阵，得到线性方程组的解x。/n/n2. 求解方程组的程序/n/nmatlab/nfunction x=Gauss(A,b)/n[n,n] = size(A); /nx = zeros(n,1); /nAug = [A,b]; %增广矩阵/n /nfor k = 1:n-1 /n    [piv,r] = max(abs(Aug(k:n,k))); %找列主元所在子矩阵的行r/n    r = r + k - 1; % 列主元所在大矩阵的行/n    if r>k /n        Aug([k,r],:)=Aug([r,k],:); /n    end/n    if Aug(k,k)==0 /n        error('对角元出现0');/n    end/n    % 把增广矩阵消元成为上三角/n    for p = k+1:n /n        Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); /n    end/nend/n /n% 解上三角方程组/nA = Aug(:,1:n); /nb = Aug(:,n+1); /nx(n) = b(n)/A(n,n); /nfor k = n-1:-1:1 /n    x(k) = (b(k)-A(k,n:-1:k+1)*x(n:-1:k+1))/A(k,k);/nend/nend/n/n% 3. 程序/nclc;/nclear;/nA   = [10 7 8 7; 7 5 6 5; 8 6 10 9; 7 5 9 10];/nb   =[32;23;33;31];/nAn =[10 7 8.1 7.2; 7.08 5.04 6 5; 8 5.98 9.89 9; 6.99 5 9 9.98];/nx    = Gauss(A,b)/nxn  = Gauss(An,b)/ndisp('The det and eigvalues of A  are:')/ndet(A),eig(A)/ndisp('The 2-norm of A  is:')/ncond(A,2)/ndisp('The value of dx  is:')/ndx =xn-x/ndisp('The 2-norm of dx  is:')/nnorm(dx,2)/ndisp('The relative error of x  is:')/nnorm(dx,2)/norm(x,2)/ndisp('The relative turbulance of A  is:')/nnorm(An-A,2)/norm(A,2)/n/n/n该代码使用Gauss函数实现高斯消元法，并计算线性方程组解的误差和扰动。/n/n* Gauss函数实现列主元高斯消元法。/n* 代码最后计算了系数矩阵A的行列式、特征值、2范数以及解的误差和扰动。/n/n### 迭代法求解线性方程组/n/n教材P211 计算实习题1/n/n参考答案：/n/n1. 雅可比迭代法求解程序/n/nmatlab/nclear;/nclc;/nn = 3;        % 可换为8和10 /nH = hilb(n)/nb = H * ones(n, 1);  /ne = 0.00001;  /nfor i = 1:n /n    if H(i, i)==0          /n        disp('对角元为零，不能求解');/n    end/nend/n%x = zeros(n, 1); /nx = ones(n,1)+rand(n,1);/nk = 0; /nkmax = 1000; /nr = 1;  /nwhile k<=kmax & r>e    /n    x0 = x;      /n    for i = 1:n         /n        s = 0;          /n        for j = 1:i - 1 /n            s = s + H(i, j) * x0(j);         /n        end/n        for j = i + 1:n /n            s = s + H(i, j) * x0(j);          /n        end/n        x(i) = b(i) / H(i, i) - s / H(i, i);     /n    end/n    r = norm(x - x0, inf); /n    k = k + 1; /nend/nif k>kmax  /n    disp('迭代不收敛，失败');     /nelse/n    disp('求解成功');/n    x     /n    k     /nend/n/n/n该代码使用雅可比迭代法求解线性方程组。/n/n* 代码首先定义系数矩阵H和常数向量b。/n* 然后初始化迭代初始值x，并设置迭代次数上限kmax和精度e。/n* 代码使用循环进行迭代计算，直到达到精度要求或达到迭代次数上限。/n/n2. SOR迭代法求解程序：/n/nmatlab/nfunction x = SOR(n, w); /nH = hilb(n);  /nb = H*ones(n, 1);  /ne = 0.00001;  /nfor i = 1:n    /n    if H(i,i)==0   /n        disp('对角线为零，不能求解');/n    end/nend/nx = zeros(n, 1); /nk = 0;/nkmax = 10000; /nr = 1; /nwhile k<=kmax & r>e  /n    x0 = x; /n    for i = 1:n     /n        s = 0;     /n        for j = 1:i - 1/n            s = s + H(i, j) * x(j);   /n        end/n        for j = i + 1:n /n            s = s + H(i, j) * x0(j);      /n        end/n        x(i) = (1 - w) * x0(i) + w / H(i, i) * (b(i) - s);/n    end/n    r = norm(x - x0, inf);   /n    k = k + 1; /nend/nif k>kmax  /n    disp('迭代不收敛，失败');    /nelse/n    disp('求解成功');/n    r/nend/nend/n/n/n该代码使用SOR迭代法求解线性方程组。/n/n* 代码首先定义系数矩阵H、常数向量b和松弛因子w。/n* 然后初始化迭代初始值x，并设置迭代次数上限kmax和精度e。/n* 代码使用循环进行迭代计算，直到达到精度要求或达到迭代次数上限。/n/n## 总结/n/n本文提供了多个数值分析代码的示例，并对每个代码段进行了详细的解释。这些代码涵盖了多项式插值、曲线拟合、数值积分以及线性方程组的直接解法和迭代解法，可以作为学习数值分析的参考。/n/n## 建议/n/n* 您可以尝试修改代码中的参数，观察结果的变化。/n* 您可以尝试使用其他数值分析方法，例如牛顿-科特斯公式、龙贝格算法、雅可比方法、高斯-赛德尔方法等。/n* 您可以尝试使用MATLAB或其他编程语言编写自己的数值分析程序。/n/n希望本文能够对您学习和理解数值分析有所帮助。