newtonofsmatlabuse.txt

the use of the matlab in newton function.

牛顿插值法的matlab实现
function a=new(c,k)
[x,y]=size(c);
b=zeros(x,2+k);
for i=1:x
for j=1:2
b(i,j)=c(i,j);
end
end
for j=3:2+k
for i=j-1:x
b(i,j)=(b(i,j-1)-b(i-1,j-1))/(b(i,1)-b(i-j+2,1));
end
end
a=b


function fx=newton(c,x,k)
b=new(c,k);
[p,q]=size(b);

fx=0; 
for i=1:p
t=b(i,i+1)
for j=1:i-1
t=t*(x-b(i-1,1));
end
fx=fx+t;
end


CCCCCCCCCCCCCCCCCCCCCCCCCCC
/****************************************************** 
* 用牛顿插值法依据N个已知数据点即使函数值 
* 输入: n--已知数据点的个数N-1 
*       x--已知数据点第一坐标的N维列向量 
* y--已知数据点第二坐标的N维列向量 
* xx-插值点第一坐标 
* 输出: 函数返回值所求插值点的第二坐标 
******************************************************/ 
double newton(int n,double x[N],double y[N],double xx) 
{ 
double d[N],b; 
int i,j; 

for(i=0;i<=n;i++) 
d=y; 
for(i=n-1;i>=0;i--) /*求差商*/ 
for(j=i+1;j<=n;j++) 
{ 
if(fabs(x-x[j])<EPS) 
return 0; 
d[j]=(d[j-1]-d[j])/(x-x[j]); 
} 
b=d[n]; 
for(i=n-1;i>=0;i--) 
b=d+(xx-x)*b; 
return b; 
} 









-----------------------------------------------
/******************************************************************** 
* 用龙贝格法计算函数f(x)从a到b的积分值 
* 输入: f--函数f(x)的指针 
*       a--积分下限 
*       b--积分上限 
*       eps--计算精度 
*       max_it--最大迭代次数 
* 输出: 返回值为f(x)从a点到b点的积分值 
*******************************************************************/ 
double romberg(double(*f)(double),double a,double b, 
                          double eps,int max_it) 
{ 
double *t,h; 
int i,m,k; 

if(!(t=(double *)malloc(max_it*sizeof(double)+1))) 
return(ERROR_CODE); 
h=b-a; 
t[1]=h*((*f)(a)+(*f)(b))/2.0; 
printf("%18.10e\n",t[1]); 
for(k=2;k<max_it+1;k++) 
{ 
double s,sm; 
h*=0.5; 
s=0.0; 
for(i=0;i<pow(2,k-2);i++) 
s+=(*f)(a+(2*i+1)*h); 
sm=t[1]; 
t[1]=0.5*t[1]+h*s; 
for(m=2;m<k+1;m++) 
{ 
s=t[m]; 
t[m]=t[m-1]+(t[m-1]-sm)/(pow(4,m-1)-1); 
if(m<k) 
sm=s; 
} 
for(m=1;m<k+1;m++) 
printf("%18.10e",t[m]); 
printf("\n"); 
if(fabs(t[k]-sm)<eps) 
{ 
sm=t[k]; 
free(t); 
return(sm); 
} 
} 
return(ERROR_CODE); 
} 








----------------------------------------------
Matlab中没有现成的拉格朗日插值函数，必须编写一个M文件实现拉格朗日插值。

设n个节点数据以数组x0, y0输入（注意Matlab的数组下标从1开始），m 个插值点以数

组x输入，输出数组y 为m 个插值。编写一个名为lagrange.m的M文件：

function y=lagrange(x0,y0,x);

n=length(x0);m=length(x);

for i=1:m

  z=x(i);

  s=0.0;

  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

 

 

 

 

 

 

拉格朗日插值的MATLAB程序

功能：给定 个不同节点，构造 的 次拉格朗日插值多项： function[c,l]=lagran(x,y)  

% 这里x为n个节点的横坐标所组成的向量，y为纵坐标所组成的向量。

% c为所得插值函数的系数组成的向量。

w=length(x);

n=w-1;

l=zeros(w,w);

for k=1:n+1

  v=1;

     for j=1:n+1

         if k~=j

             v=conv(v,poly(x(j)))/(x(k)-x(j));

    end

end

  l(k,:)=v;

end

c=y*l;

牛顿插值的MATLAB程序

程序功能：给定 个不同节，构造次数小于等于 的牛顿插值多项式(12.3.7)。

 

******************************************************************

function[c,d]=newpoly(x,y)  

% 这里 x为n个节点的横坐标所组成的向量，y为纵坐标所组成的向量。

% c为所求的牛顿插值多项式的系数构成的向量。

n=length(x);

d=zeros(n,n);

d(:,1)=y';

for j=2:n

    for k=j:n

     d(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));

  end

end

c=d(n,n);

for k=(n-1):-1:1

   c=conv(c,poly(x(k)));

   m=length(c);

   c(m)=c(m)+d(k,k);

end