插值法,迭代法matlab程序

合集下载
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

数值分析作业

姓名王建忠

学号132080202006

学院能源与动力工程

专业机械电子工程

2013年12月16日

1.用MATLAB编程实现langrage插值多项式:

syms x

x0=[-2,-1,0,1];

y0=[3,1,1,6];

n=length(x0);

for i=1:n

a=1;

for j=1:n

if j~=i

a=expand(a*(x-x0(j)));

end

end

b=1;

for k=1:n

if k~=i

b=b*(x0(i)-x0(k));

end

end

A(i)=expand(a/b);

end

L=0;

for p=1:n

L=L+y0(p)*A(p);

end

L

>> Language

L =x^3/2 + (5*x^2)/2 + 2*x + 1

2.牛顿插值多项式程序

function [p2,z]=newTon(x,y,t)

%输入参数中x,y为元素个数相等的向量,t为待估计的点,可以为数字或向量。%输出参数中p2为所求得的牛顿插值多项式,z为利用多项式所得的t的函数值。

n=length(x);

chaS(1)=y(1);

for i=2:n

x1=x;y1=y;

x1(i+1:n)=[];

y1(i+1:n)=[];

n1=length(x1);

s1=0;

for j=1:n1

t1=1;

for k=1:n1

if k==j

continue;

else

t1=t1*(x1(j)-x1(k));

end

end

s1=s1+y1(j)/t1;

end

chaS(i)=s1;

end

b(1,:)=[zeros(1,n-1) chaS(1)];

cl=cell(1,n-1);

for i=2:n

u1=1;

for j=1:i-1

u1=conv(u1,[1 -x(j)]);

cl{i-1}=u1;

end

cl{i-1}=chaS(i)*cl{i-1};

b(i,:)=[zeros(1,n-i),cl{i-1}];

end

p2=b(1,:);

for j=2:n

p2=p2+b(j,:);

end

if length(t)==1

rm=0;

for i=1:n

rm=rm+p2(i)*t^(n-i);

end

z=rm;

else

k1=length(t);

rm=zeros(1,k1);

for j=1:k1

for i=1:n

rm(j)=rm(j)+p2(i)*t(j)^(n-i);

end

z=rm;

end

end

x=[-1 1 2 5];

y=[-7 7 -4 35];

t=1;

[u,v]=newTon(x,y,t)

>> [u,v]=newton(x,y,t)

u =2 -10 5 10

v=7

则所求得的牛顿多项式为:2*x^3-10x^2+ 5*x +10

3.Matlab程序三次样条插值函数

clear

clc

x=[1 2 3 4 6]

y=[0 0.693147 1.09861 1.38629 1.79176]

n=length(x);

for i=1:n-1

h(i)=x(i+1)-x(i);

end

for i=1:n-2

k(i)=h(i+1)/(h(i)+h(i+1));

u(i)=h(i)/(h(i)+h(i+1));

end

for i=1:n-2

gl(i)=3*(u(i)*(y(i+2)-y(i+1))/h(i+1)+k(i)*(y(i+1)-y(i))/h(i)); end

g0=3*(y(2)-y(1))/h(1);

g00=3*(y(n)-y(n-1))/h(n-1);

g=[g0 gl g00];

g=transpose(g)

k1=[k 1];

u1=[1 u];

Q=2*eye(5)+diag(u1,1)+diag(k1,-1)

m=transpose(Q\g)

syms X;

for i=1:n-1

p1(i)=(1+2*(X-x(i))/h(i))*((X-x(i+1))/h(i))^2*y(i);

p2(i)=(1-2*(X-x(i+1))/h(i))*((X-x(i))/h(i))^2*y(i+1);

p3(i)=(X-x(i))*((X-x(i+1))/h(i))^2*m(i);

p4(i)=(X-x(i+1))*((X-x(i))/h(i))^2*m(i+1);

p(i)=p1(i)+p2(i)+p3(i)+p4(i);

p(i)=simple(p(i));

end

s1=p(1)

相关文档
最新文档