常用数值算法matlab源码(pde见另贴)() 下载本文

%举例:

%x=[0.32 0.34 0.36];y=[0.314567 0.333487 0.352274];x0=0.33; %f=InterpLagrange(x,y) %f=InterpLagrange(x,y,x0)

if length(x)==length(y) n=length(x); else

disp('节点个数和函数值个数不同!') f=' '; return; end p=0; for i=1:n l=y(i); for j=1:n if j==i continue; end

%利用卷积计算Lagrange基函数 l=conv(l,[1 -x(j)]./(x(i)-x(j))); end

%p是一向量,表示插值多项式的系数 p=p+l;

end

if nargin==3

f=polyval(p,x0);%计算插值多项式在x0处的值 else

f=poly2str(p,'x');%把插值多项式的向量形式转化为插值多项式的符号形式 end

function f=InterpLagrange2(x,y,x0) %构造Lagrange插值多项式

%此函数中借助符号运算来求Lagrange基函数,运算速度较慢,不推荐此种方法 %f=InterpLagrange2(x,y,x0)

%f:插值多项式或者是插值多项式在x0处的值 %x:节点 %y:函数值 %x0:某一测试点 %

%调用格式:

%f=InterpLagrange2(x,y) 返回插值多项式

%f=InterpLagrange2(x,y,x0) 返回插值多项式在点x0处的值 %举例:

%x=[0.32 0.34 0.36];y=[0.314567 0.333487 0.352274];x0=0.33; %f=InterpLagrange2(x,y)

%f=InterpLagrange2(x,y,x0)

if length(x)==length(y) n=length(x); else

disp('节点个数和函数值个数不同!') f=' '; return; end syms t; f=0; for i=1:n l=y(i); for j=1:n if j==i continue; end

l=l*(t-x(j))/(x(i)-x(j));%借助符号运算,计算Lagrange基函数 end f=f+l;

simplify(f);%化简多项式 if i==n

if nargin==3

f=subs(f,'t',x0);%计算插值多项式f在点x0处的值 else

f=collect(f);%计算插值多项式,展开并合并同类项 f=vpa(f,6);%设置多项式系数的有效数字 end end end

4、Newdon插值法

function f=InterpNewdon(x,y,x0) %Newdon插值多项式 %f=InterpNewdon(x,y,x0)

%f:插值多项式或者是插值多项式在x0处的值 %x:节点 %y:函数值 %x0:某一测试点 %

%调用格式

%f=InterpNewdon(x,y) 返回插值多项式

%f=InterpNewdon(x,y,x0) 返回插值多项式在x0点的值 %应用举例:

%x=[1 2 3 4 5];y=[1 4 7 8 6];x0=6;

%f=InterpNewdon(x,y) %f=InterpNewdon(x,y,x0)

if length(x)==length(y) n=length(x); else

disp('节点个数和函数值个数不同!') f=' '; return; end

A=zeros(n);%初始化差商矩阵 for i=1:n

A(i,1)=y(i);%差商矩阵的第一列是函数值 end

%计算差商矩阵

%差商矩阵中对角线上的元素为Newdon插值多项式的系数 for j=2:n for i=j:n

A(i,j)=(A(i,j-1)-A(i-1,j-1))/(x(i)-x(i-j+1)); end end

%求Newdon插值多项式 p=zeros(1,n);