数值分析实验题和程序

一、实验3.1

编制以函数{}

n k

k x

=危机的多项式最小二乘拟合程序，并用于对表3.11中的数据作3次

取权数1i w ≡，求拟合曲线3

*

k

k k a x

φ==

∑中的参数{}k a 、平方误差2

σ，并作离散数据

{},i i x y 的拟合曲线*()y x φ=的图形。

程序代码： x0=-1:0.5:2;

y0=[-4.447 -0.452 0.551 0.048 -0.447 0.549 4.552]; n=3;

alph=polyfit(x0,y0,n) %参数{a k } y=polyval(alph,x0);

r=(y0-y)*(y0-y)' %平方误差2

σ=r y=polyval(alph,x); x=-1:0.01:2; y=plot(x,y,'k--'); xlabel('x');

ylabel('y0 * and polyfit.y--'); hold on;

plot(x0,y0,'*');

title('离散数据的3项拟合') gridon;

实验结果：

拟合函数*

()y x φ=的图形：

拟合曲线

3

*

k

k

k

a x

φ

=

=∑中的参数{}k a中3210

,,,

a a a a依次为alph中的四个数值。

alph =

1.9991 -

2.9977 -0.0000 0.5491

平方误差2

σ=r。

r =

2.1762e-005

实验分析：

最小二乘曲线拟合是在离散情形下的最佳平方逼近，拟合的曲线很光滑，而且所有的7个数值点均在曲线上，拟合效果很好；拟合的平方误差很小，为10-5量级。

二、实验3.2

编制正交化多项式最小二乘拟合程序，并用于求解上题中的3次多项式最小二乘拟合问题，作拟合曲线的图形，计算平方误差，并与上题结果进行比较。

程序代码：

x=-1:0.5:2;

y=[-4.447 -0.452 0.551 0.048 -0.447 0.549 4.552];

n=3;

result=inputdlg({'请输入权向量w:'},'charpt-3',1,{'[1 1 1 1 1 1 1]'});

w=str2num(char(result));

m=length(x)-1;

s1=0;s2=ones(1,m+1);v2=sum(w);

d(1)=y*w';c(1)=d(1)/v2;

for k=1:n

xs=x.*s2.^2*w';a(k)=xs/v2;

if(k==1)

b(k)=0;

else

b(k)=v2/v1;

end

s3=(x-a(k)).*s2-b(k)*s1;

v3=s3.^2*w';

d(k+1)=y.*s3*w';c(k+1)=d(k+1)/v3;

s1=s2;s2=s3;v1=v2;v2=v3;

end

r=y.*y*w'-c*d'

alph=zeros(1,n+1)

T=zeros(n+1,n+2);

T(:,2)=ones(n+1,1);T(2,3)=-a(1);

if(n>=2)

for k=3:n+1

for i=3:k+1

T(k,i)=T(k-1,i)-a(k-1)*T(k-1,i-1)-b(k-1)*T(k-2,i-2);

end

end

end

for i=1:n+1

for k=i:n+1

alph(n+2-i)=alph(n+2-i)+c(k)*T(k,k+2-i);

end

end

xmin=min(x);xmax=max(x);dx=(xmax-xmin)/(25*m);

t=(xmin-dx):dx:(xmax+dx);

s=alph(1);

for k=2:n+1

s=s.*t+alph(k);

end

plot(x,y,'x',t,s,'-');

xlabel('x');

ylabel('y');

grid on;

disp(alph);

disp(r);

实验结果：

拟合曲线图形：

参数{}k a 中3210,,,a a a a 依次为alph 中的四个数值：1.9991 -2.9977 -0.0000 0.5491 平方误差2

σ=r ：2.1762e-005

实验分析：

比较实验3.1和3.2的结果发现：对于同一个数据表，两种方法的拟合参数、误差均是相等的，表示这两种方法的拟合效果是一样的。因为这个数据表不是病态的法方程组，所以没有体现出正交多项式曲线拟合的优点。根据这两个实验代码和运算时间来看，我们在做曲线拟合时可以先初步判断数据点是否病态，若否，就直接用最小二乘法拟合。以避免不必要的程序复杂化。

三、实验5.1 常微分方程性态和R-K 法稳定性试验

实验目的：考察下面微分方程右端项中函数y 前面的参数对方程性态的影响（它可以使方程为好条件的或坏条件的）和研究计算步长对R-K 法计算稳定性的影响。 实验题目：常微分方程初值问题

'1,01,(0)1,

y ay ax x y ⎧=-+<<⎨

=⎩其中，5050a -≤≤.其精确解为()ax

y x e x =+。 实验要求：本实验题都用4阶经典R-K 法计算。

（1）对参数a 分别取4个不同的数值：一个大的正值，一个小的正值，一个绝对值小的负值和一个绝对值大的负值。取步长h=0.01，分别用经典的R-K 法计算，将四组计算结果画在同一张图上，进行比较并说明相应初值的性态。

（2）取参数a 为一个绝对值不大的负值和两个计算步长，一个步长使参数ah 在经典R-K 法的稳定域内，另一个步长在经典R-K 法的稳定域外。分别用经典R-K 法计算并比较结果。取全域等距的10个点上的计算值，列表说明。

程序代码（1）： function charp5_1 clf;