插值法,迭代法matlab程序

数值分析作业

姓名王建忠

学号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)

s2=p(2)

s3=p(3)

s4=p(4)

for k=1:4

for z=x(k):0.001:x(k+1)

q=eval(subs(p(k),'X','z'));

plot(z,q,'b')

hold on

end

end

grid on

legend('eval')

title('sysm')

xlabel('x')

ylabel('p')

x =

1 2 3 4 6

y =

0 0.6931 1.0986 1.3863 1.7918

g = 2.0794

1.6479

1.0397

0.7781

0.6082

Q = 2.0000 1.0000 0 0 0

0.5000 2.0000 0.5000 0 0

0 0.5000 2.0000 0.5000 0

0 0 0.6667 2.0000 0.3333

0 0 0 1.0000 2.0000

m = 0.7629 0.5537 0.3181 0.2535 0.1774

4.牛顿迭代法求多项式根

syms x

f=x-cos(x);

df=diff(f,x);

eps=1e-6;

x0=1;

cnt=0;

MAXCNT=20; %最大循环次数

while cnt

x1=x0-subs(f,x,x0)/subs(df,x,x0); %去掉这个分号,可以看到迭代过程.

if (abs(x1-x0)

break;

end

x0=x1;

cnt=cnt+1;

end

if cnt==MAXCNT

disp '不收敛'

else

vpa(x1,8)

end

>> Untitled

ans =

0.73908513