Matlab数值计算示例： 牛顿插值法、LU分解法、拉格朗日插值法、牛顿插值法

本文源于一次课题作业，部分自己写的，部分借用了网上的demo

• 牛顿迭代法（1）

x=1:0.01:2;y=x.^3-x.^2+sin(x)-1;plot(x,y,'linewidth',2);grid on;%由图像可知 根在1.05到1.15之间syms xs0=diff(x^3-x^2+sin(x)-1,x,1);% 得到s0= cos(x) - 2*x + 3*x^2% 迭代方程为 y=x-(x.^3-x.^2+sin(x)-1)/(cos(x) - 2.*x + 3*x.^2)clear;x=1.15;for i=1:30    x=x-(x.^3-x.^2+sin(x)-1)/(cos(x) - 2.*x + 3*x.^2)%根据牛顿迭代法公式。一直迭代计算30次。end%可得  x取值为1.0935;

• 牛顿迭代法（2）

%%   绘制图形。判断跟的大概位置。x=1:0.01:2;f=x.^3-x.^2+sin(x)-1;plot(x,f,'linewidth',2);grid on;%由图像可知 根在1.05到1.15之间%%clc,clearsyms x f=x.^3-x.^2+sin(x)-1;%所求函数df=diff(f,x); %求取一阶导数eps=1e-4; %误差判断x0=1.15; %迭代初始值。cnt=0; MAXCNT=20; %最大循环次数 while cnt<MAXCNT %防止无限循环        x1=x0-subs(f,x,x0)/subs(df,x,x0); %去掉这个分号,可以看到迭代过程.      if (abs(x1-x0)<eps)       break;       end         x0=x1;         cnt=cnt+1; end if cnt==MAXCNT     disp '不收敛' else     vpa(x1,8) end

• LU分解法

被调函数:

function [L,U]=lufj(A)%  利用紧凑格式法原理 编写的LU 分解[n,m]=size(A); % 获取A矩阵的行和列if m~=n    %判断行列相等与否    error('Not a squared matrix1');else    A(2:n)=A(2:n)/A(1,1);    for k=2:n-1        A(k,k:n)=A(k,k:n)-A(k,1:k-1)*A(1:k-1,k:n);        A(k+1:n,k)=(A(k+1:n,k)-A(k+1:n,1:k-1)*A(1:k-1,k))/A(k,k);    end    %都是根据定义进行循环计算。endL=A;U=A;for i=1:n    L(i,i)=1;    L(i,i+1:n)=0;    U(i,1:i-1)=0;end

主函数：

%%  需要调用lufj函数；A=[-2 -2 3 5    1 2 1 -2    2 5 3 -2    1 3 2 3]b=[-1    4    7    0]% x=A\b  %左除法求解[L,U]=lufj(A);x0=L\b;x=U\x0%求出的x即为解   

• 拉格朗日插值法
  被调函数：

function y=lagrange(x0,y0,x);% 根据拉格朗日插值定义编写n=length(x0);m=length(x);for i=1:m  z=x(i);  s=0.0;%给s的初值  for k=1:n      p=1.0;      for j=1:n          if j~=k            p=p*(z-x0(j))/(x0(k)-x0(j));          end      end      s=p*y0(k)+s;    end    y(i)=s;end

主函数：

x=[0,1,2,4];y=[1,9,23,3];y0=lagrange(x,y,1.5)

• 牛顿插值
  被调函数：

function yi=New_Int(x,y,xi)%Newton基本插值公式%x为向量，全部的插值节点%y为向量，差值节点处的函数值%xi为标量，是自变量%yi为xi出的函数估计值n=length(x);m=length(y);if n~=m    error('The lengths of X ang Y must be equal!');    return;end%计算均差表YY=zeros(n);Y(:,1)=y';for k=1:n-1    for i=1:n-k        if abs(x(i+k)-x(i))<eps            error('the DATA is error!');            return;        end        Y(i,k+1)=(Y(i+1,k)-Y(i,k))/(x(i+k)-x(i));    endend%计算牛顿插值公式yi=0;for i=1:n    z=1;    for k=1:i-1        z=z*(xi-x(k));    end    yi=yi+Y(1,i)*z;end

主函数：

clear all clcx0=[0.4 0.55 0.65 0.80 0.90 1.05];y0=[0.41075 0.57815 0.69675 0.88811 1.0265 1.25382]; x1=0.596;  % 待插值点。y1=New_Int(x0,y0,x1)% y1即为待插值点的函数值。

TIP：主函数和被调函数要放在一个文件夹内。否则会引起调用错误

NOTE:本文对基本方法做了总结，你可以结合理论知识再来看代码，希望对你有所帮助