数值分析(第五版)计算实习题第三章

数值分析计算实习题第三章

第二次作业：

题一：

x=-1:0.2:1;y=1./(1+25.*x.^2);

f1=polyfit(x,y,3)

f=poly2sym(f1)

y1=polyval(f1,x)

x2=linspace(-1,1,10)

y2=interp1(x,y,x2)

plot(x,y,'r*-',x,y1,'b-')

hold on

plot(x2,y2,'k')

legend('数据点','3次拟合曲线','3次多项式插值')

xlabel('X'),ylabel('Y')

输出：f1 =

0.0000 -0.5752 0.0000 0.4841

f =

(4591875547102675*x^3)/81129638414606681695789005144064 - (3305*x^2)/5746 + (1469057404776431*x)/20282409603651670423947251286016 + 4360609662300613/9007199254740992

y1 =

-0.0911 0.1160 0.2771 0.3921 0.4611 0.4841 0.4611 0.3921 0.2771 0.1160 -0.0911

x2 =

-1.0000 -0.7778 -0.5556 -0.3333 -0.1111 0.1111 0.3333 0.5556 0.7778 1.0000

y2 =

0.0385 0.0634 0.1222 0.3000 0.7222 0.7222 0.3000 0.1222 0.0634 0.0385

题二：

X=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];

Y=[1.0 0.41 0.50 0.61 0.91 2.02 2.46];

p1=polyfit(X,Y,3)

p2=polyfit(X,Y,4)

Y1=polyval(p1,X)

Y2=polyval(p2,X)

plot(X,Y,'r*',X,Y1,'b-.',X,Y2,'g--')

p3=polyfit(X,Y,2)

Y3=polyval(p3,X)

f1=poly2sym(p1)

f2=poly2sym(p2)

f3=poly2sym(p3)

plot(X,Y,'r*',X,Y1,'b-.',X,Y2,'g--',X,Y3,'m--')

legend('数据点','3次多项式拟合','4次多项式拟合','2次多项式拟合') xlabel('X轴'),ylabel('Y轴')

输出：

p1 =

-6.6221 12.8147 -4.6591 0.9266

p2 =

2.8853 -12.3348 16.2747 -5.2987 0.9427

Y1 =

0.9266 0.5822 0.4544 0.5034 0.9730 2.0103 2.4602 Y2 =

0.9427 0.5635 0.4399 0.5082 1.0005 1.9860 2.4692 p3 =

3.1316 -1.2400 0.7356

Y3 =

0.7356 0.6429 0.6128 0.6454 0.8984 1.7477 2.6271

f1 =

- (7455778416425075*x^3)/1125899906842624 + (1803512222945435*x^2)/140737488355328 - (40981580032809*x)/8796093022208 + 8345953784399011/9007199254740992

f2 =

(1624271450198125*x^4)/562949953421312 - (3471944732519173*x^3)/281474976710656 + (4580931990070659*x^2)/281474976710656 - (1491459232922115*x)/281474976710656 + 1061409433081293/1125899906842624

f3 =

(18733*x^2)/5982 - (74179*x)/59820 + 73337/99700

题三：

建立三角插值函数的m文件

function [A,B,Y1,Rm]=sanjiaobijin(X,Y,X1,m)%A B分别是m阶三角多项式Tm（x）的系数aj，bj（j=1,2,...,m）的系数矩阵，Y1是Tm（x）在X1处的值，X Y数据点,Rm为均方误差

n=length(X)-1;max1=fix((n-1)/2);

if m>max1

m=max1;

end

A=zeros(1,m+1);B=zeros(1,m+1);

Ym=(Y(1)+Y(n+1))/2;Y(1)=Ym;

Y(n+1)=Ym;A(1)=2*sum(Y)/n;

for i=1:m

B(i+1)=sin(i*X)*Y';

A(i+1)=cos(i*X)*Y';

end

A=2*A/n;B=2*B/n;

A(1)=A(1)/2;Y1=A(1);

for k=1:m

Y1=Y1+A(k+1)*cos(k*X1)+B(k+1)*sin(k*X1);

Tm=A(1)+A(k+1).*cos(k*X)+B(k+1).*sin(k*X);k=k+1;

end

Y,Tm,Rm=(sum(Y-Tm).^2)/n

输出：>> X=-pi:2*pi/33:pi;