数值分析实验题和程序

一、实验3.1编制以函数{}n kk x=危机的多项式最小二乘拟合程序，并用于对表3.11中的数据作3次取权数1i w ≡，求拟合曲线3*kk 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*kkka 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:nxs=x.*s2.^2*w';a(k)=xs/v2;if(k==1)b(k)=0;elseb(k)=v2/v1;ends3=(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;endr=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+1for i=3:k+1T(k,i)=T(k-1,i)-a(k-1)*T(k-1,i-1)-b(k-1)*T(k-2,i-2);endendendfor i=1:n+1for k=i:n+1alph(n+2-i)=alph(n+2-i)+c(k)*T(k,k+2-i);endendxmin=min(x);xmax=max(x);dx=(xmax-xmin)/(25*m);t=(xmin-dx):dx:(xmax+dx);s=alph(1);for k=2:n+1s=s.*t+alph(k);endplot(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 -≤≤.其精确解为()axy 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;result=inputdlg({'请输入参数a:'},'charpt-5',1,{'40'});a=str2num(char(result));result=inputdlg({'请输入（0 1）之间的步长:'},'实验5.1',1,{'0.01'}); h=str2num(char(result));x=0:h:1;y=x;N=length(x);y(1)=1;func=inline('1+(y-x).*a');for n=1:N-1k1=func(a,x(n),y(n));k2=func(a,x(n)+h/2,y(n)+k1*h/2);k3=func(a,x(n)+h/2,y(n)+k2*h/2);k4=func(a,x(n)+h,y(n)+k3*h);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endy0=log(exp(a*x)+x);result=inputdlg({'请输入参数b:'},'charpt-5',1,{'10'});b=str2num(char(result));result=inputdlg({'请输入（0 1）之间的步长:'},'实验5.1',1,{'0.01'}); h=str2num(char(result));x=0:h:1;y=x;N=length(x);y(1)=1;func=inline('1+(y-x).*b');for n=1:N-1k1=func(b,x(n),y(n));k2=func(b,x(n)+h/2,y(n)+k1*h/2);k3=func(b,x(n)+h/2,y(n)+k2*h/2);k4=func(b,x(n)+h,y(n)+k3*h);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endy1=log(exp(b*x)+x);result=inputdlg({'请输入参数c:'},'charpt-5',1,{'-10'});c=str2num(char(result));result=inputdlg({'请输入（0 1）之间的步长:'},'实验5.1',1,{'0.01'}); h=str2num(char(result));x=0:h:1;y=x;N=length(x);y(1)=1;func=inline('1+(y-x).*c');for n=1:N-1k1=func(c,x(n),y(n));k2=func(c,x(n)+h/2,y(n)+k1*h/2);k3=func(c,x(n)+h/2,y(n)+k2*h/2);k4=func(c,x(n)+h,y(n)+k3*h);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endy2=log(exp(c*x)+x);result=inputdlg({'请输入参数d:'},'charpt-5',1,{'-40'});d=str2num(char(result));result=inputdlg({'请输入（0 1）之间的步长:'},'实验5.1',1,{'0.01'}); h=str2num(char(result));x=0:h:1;y=x;N=length(x);y(1)=1;func=inline('1+(y-x).*d');for n=1:N-1k1=func(d,x(n),y(n));k2=func(d,x(n)+h/2,y(n)+k1*h/2);k3=func(d,x(n)+h/2,y(n)+k2*h/2);k4=func(d,x(n)+h,y(n)+k3*h);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endy3=log(exp(d*x)+x);plot(x,y0,x,y1,x,y2,x,y3);hold on;xlabel('x');ylabel('y=log(exp(a*x)+x)');legend('y0(a=40)','y1(a=10)','y2(a=-10)','y3(a=-40)');实验结果（1）：注释：由于y 轴数据差异较大，所以图形中对y 值取对数。

实验分析（1）：经典R-K 算法的绝对稳定区域是 2.785h 0a -≤≤，因此当h=0.01时，a 的取值范围应该是278.50a -≤≤ 但有5050a -≤≤，因此a 的稳定区域为500a -≤≤。