计算方法实验指导书

1、计算方法实验指导书数学实验室2013年3月目 录实验一病态问题2实验二误差传播与算法稳定性6实验三线方程组的直接解法9实验四解线性方程组的迭代法12实验五非线性方程求根13实验六函数插值方法14提高：多项式插值的振荡现象15扩展：样条插值的收敛19实验七数值积分与数值微分23实验八常微分方程初值问题数值解法24实验九函数逼近与曲线拟合26附实验一：多项式最小二乘拟合27附实验二：正交化多项式的最小二乘拟合29附录一：实验报告内容要求32附录二：部分程序示例331、线方程组的直接解法系列程序332、Jacobi迭代法解线性方程组程序383、牛顿法非线性方程求根程序394、多项式系列程序405、R

2、omberg算法计算数值积分程序426、四阶经典的龙格-库塔方法解常微分方程程序427、最小二乘法曲线拟合程序438、幂法求矩阵特征值程序46前 言结合课程教学，配备适当的上机实验(16个学时)以便加深课堂教学的实践性，同时通过实验可以加强对数学模型的总体分析，算法选取，程序结构，上机调试和结果分析等环节的训练。本实验指导书共包含8个实验，要求学生在16个实验课时内完成。为使实验更为有成效，需要写出实验报告（格式要求见附录）,以此可作为计算方法课程成绩评定的参考。实验一病态问题一、实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题是指问题本身对扰动

3、敏感，反之属于好问题。本实验通过对一个高次多项式方程的求解，初步认识病态问题。二、实验设备和实验环境安装有C、C+或MATLAB的计算机。三、实验内容及要求：考虑一个高次的代数多项式 (1.1)显然该多项式的全部根为1，2，20，共计20个，且每个根都是单重的(也称为简单的)。现考虑该多项式的一个扰动 (1.2)其中，是一个非常小的数。这相当于是对方程(1)中x19的系数作一个小的扰动。比较方程(1.1)和方程(1.2)根的差别，从而分析方程(1.1)的解对扰动的敏感性。四、实验步骤与结果分析(一) 实验源程序function t_charpt1_1% 数值实验1病态问题% 输入：0 20之间

4、的扰动项及小的扰动常数% 输出：加扰动后得到的全部根clcresult=inputdlg('请输入扰动项:在0 20之间的整数:','charpt 1_1',1,'19');Numb=str2num(char(result);if(Numb>20)|(Numb<0)errordlg('请输入正确的扰动项:0 20之间的整数!');return;endresult=inputdlg('请输入(0 1)之间的扰动常数:','charpt 1_1',1,'0.00001');e

5、ss=str2num(char(result);ve=zeros(1,21);ve(21-Numb)=ess;root=roots(poly(1:20)+ve);x0=real(root); y0=imag(root);plot(x0',y0', '*');disp('对扰动项 ',num2str(Numb),'加扰动',num2str(ess),'得到的全部根为:');disp(num2str(root);(二) 实验结果分析(1) 对于x19项的扰动ess，不同的取值对应的结果如下所示。l 对扰动项 19加扰动

6、1e-010得到的全部根为:19.9961，19.0257，17.9085，17.1508，15.7982，15.181，13.8995，13.0571，11.9753，11.0109，9.99608，9.00111，7.99978，7.00003，6，5，4，3，2，1。l 对扰动项 19加扰动1e-009得到的全部根为:19.952，19.2293，17.6573+0.692896i，17.6573-0.692896i，15.4524+0.875524i，15.4524-0.875524i，13.3527+0.486992i，13.3527-0.486992i，11.8578，11.042

7、7，9.9916，9.00201，7.99952，7.00009，5.99999，5，4，3，2，1。l 对扰动项 19加扰动1e-007得到的全部根为:20.422+0.999203i，20.422-0.999203i，18.1572+2.4702i，18.1572-2.4702i，15.3149+2.69865i，15.3149-2.69865i，12.8466+2.06246i，12.8466-2.06246i，10.9216+1.10366i，10.9216-1.10366i，9.56629，9.11508，7.99387，7.00027，6，5，4，3，2，1。l 对扰动项 19加扰

8、动1e-005得到的全部根为:22.5961+2.3083i，22.5961-2.3083i，18.8972+5.00563i，18.8972-5.00563i，14.9123+4.95848i，14.9123-4.95848i，12.0289+3.73551i，12.0289-3.73551i，10.059+2.33021i，10.059-2.33021i，8.63828+1.0564i，8.63828-1.0564i，7.70896，7.028，5.99942，5.00001，4，3，2，1。根在复平面上的位置如图所示：图 ess=1e-010 图 ess=1e-009图 ess= 1e-

9、007 图 ess=1e-005从实验的图形中可以看出，当ess充分小时，方程1.1和方程1.2的解相差很小，当ess逐渐增大时，方程的解就出现了病态解，这些解都呈现复共轭性质。并且，病态解首先出现在x=16这个解附近，如ess=1e-009时，x=20，19，12，11，2，1的解基本误差不大。在x=16附近，扰动后的解偏离实轴程度较严重，随着ess的增大，扰动对解的影响从x=16附近开始向两边波及，并且偏离实轴的幅度越来越大。x=0，1，2，3，4，5这些阶次较小的解对x19上的扰动最不敏感。(2) 将扰动项加到x18上后，ess=1e-009时方程的解都比较准确，没有出现复共轭现象。es

10、s=1e-008时误差与x19(ess=1e-009)时相当，即扰动加到x18上比加到x19小一个数量级。对x8的扰动ess=1000时没有出现复共轭，误差很小；对x的扰动ess=10e10时没有出现复共轭，误差很小。因此，扰动作用到xn上时，n越小，扰动引起的误差越小。(3)令： (1.3)的零点均为的函数，分别它们记为， 显然有。研究关于的变化情况，将表示为： (1.4)关于的变化或敏感程度可以用其导数表示，故在（1.4）两边关于求导：(1.5)为了知道原方程的根是如何受扰动的影响，需要知道。在（1.5）两边令，得到 (1.6)再令，则得，即(1.7)由于，故(1.8)计算表明，对根1，此

11、导数的绝对值只有，极其微小；但从第7个根起，此导数的绝对值就从开始，最大直到，非常大！所以必定造成病态。这就是根源。现在来估计扰动对根的影响。对根 的影响，由（1.7）可得条件数：(1.9)经过计算发现，扰动对的影响最大，对的影响最小，与实际计算结果一致。实验二误差传播与算法稳定性一、实验目的体会稳定性在选择算法中的地位。误差扩张的算法是不稳定的，是我们所不期望的；误差衰竭的算法是稳定的，是我们努力寻求的。二、实验设备和实验环境安装有C、C+或MATLAB的计算机。三、实验内容及要求：考虑一个由积分定义的序列(2.1)显然En>0，n=1，2，当n=1时，利用部分积分得E1=1/e。而对

12、n2，经分部积分可得递推关系En=1-nEn-1, n=2,3, (2.2)由式(2.1)得En1/(n+1)。由递推关系式(2.2)，可得计算式(2.1)积分序列En的两种算法。其一为式(2.2)的直接应用，即：E1=1/e, En= 1nEn-1, n=2, 3, (2.3)另一种算法则是利用式(2.2)变形得到：EN=0，En-1=(1-En)/n, n=N-1, N-2, , 3, 2. (2.4)四、实验步骤与结果分析(一) 实验源程序function t_charpt1_2% 数值实验1.2:误差传播与算法稳定性% 输入:递推式选择及递推步数% 输出：各步递推值及误差结果，以及递推

13、值和误差与递推步数的关系图clcpromps='请选择递推关系式，若选公式2.3,请输入1，否则输入2:'result=inputdlg(promps,'charpt 1_2',1,'1');Nb = str2num(char(result);if(Nb=1)&(Nb=2) errordlg('请选择递推关系式，若选公式2.3,请输入1，否则输入2!');return;endresult=inputdlg('请输入递推步数n:','charpt 1_2',1,'10');st

14、eps=str2num(char(result);if(steps<1)errordlg('递推步数错误!');return;endresult=zeros(1,steps); err=result;if(Nb=1) n=1; result(n)=1/exp(1); while(n<steps) result(n+1)=1-n*result(n); err(n+1)=abs(result(n+1)-func(n+1); n=n+1; endelseif(Nb=2) n=steps; err(n)=abs(result(n)-func(n); while(n>1

15、) result(n-1)=(1-result(n)/n; err(n-1)=abs(result(n-1)-func(n-1); n=n-1; endenddisp('递推值:',num2str(result);disp('误差:',num2str(err);plot(1:steps,result,'-');grid onhold on;plot(1:steps,err,'r-');xlabel('n'); ylabel('En- and ERRn-');text(2,err(2),'up

16、arrow err(n)');text(4,result(4),'downarrow En');%-function en=func(n)% 计算En的精确值if(n=1) en=1/exp(1);else en=1-n*func(n-1);end(二) 实验结果分析(1)分别用算法(2.3)、算法(2.4)计算得到的结果如图所示。图 算法(2.3)结果 图 算法(2.4)结果两种算法的得到的结果数据如下：l 算法(2.3)：递推值: 0.367879441 0.632120559 -0.264241118 1.79272335 -6.1708934131.854467

17、1 -190.126802 1331.88762 -10654.1009 95887.9084误差: 0 0.367879441 0.471517765 1.62182994 6.31642635 31.7276647 190.239186 1331.78668 10654.1925 95887.8245l 算法(2.4)：递推值:0.36788 0.26424 0.20728 0.17089 0.14554 0.12679 0.1125 0.1 0.1 0误差:2.3114e-008 4.6229e-008 1.3869e-007 5.5474e-007 2.7737e-006 1.6642

18、e-0050.0001165 0.00093197 0.0083877 0.083877显然可以看出，算法(2.4)可以给出精确的结果，算法(2.4)虽然迭代的第一次EN的误差可能比较大，但EN-1的误差就急剧缩小，后面的误差都很小；而算法(2.3)的误差却是逐渐增大的，E10的误差居然达到了惊人的9.6´104。(2)两种算法的优劣与第一感觉是吻合的。理论上分析如下：算法(2.3)：En= 1nEn-1，因此误差en=|En-En*|= n|En-1-En-1*|= n! |E1-E1*|= n!e1。可见误差的传递呈n的阶乘增大，算法是不稳定的。算法(2.4)：En-1=(1-E

19、n)/n，因此误差为： 。因此，算法的误差是逐级递减的，当N取的比较大时，即使初始EN的误差很大，计算的En的误差也会很小，因此算法是稳定的。(3)算法(2.3)中即使e1很小，由于n增大时，误差呈n倍增大，因此算法(2.3)的误差对后面的影响是扩张的。算法(2.4)中即使EN很大，由于误差呈1/n递减，因此算法(2.4)的误差对后面的影响是衰减的。(4)综合上面的分析，算法(2.3)是不稳定的，算法(2.4)是稳定的。实验三线方程组的直接解法一、实验目的1掌握求解线性方程组的高斯消去法¾列选主元在计算机上的算法实现。2程序具有一定的通用性，程序运行时先输入一个数n表示方程含有的未知

20、数个数，然后输入每个线性方程的系数和常数，求出线性方程组的解。二、实验设备和实验环境安装有C、C+或MATLAB的计算机。三、实验内容：利用列选主元高斯消去法求解线性方程组：四、算法描述：1高斯消去法基本思路设有方程组Ax = b，设A是可逆矩阵。高斯消去法的基本思想就是将矩阵的初等行变换作用于方程组的增广矩阵B = A»b，将其中的A变换成一个上三角矩阵，然后求解这个三角形方程组。2列主元高斯消去法计算步骤将方程组用增广矩阵表示。步骤1：消元过程，对k = 1, 2, , n-1(1) 选主元，找ikÎk, k+1, , n使得：(2) 如果，则矩阵A奇异，程序结束；否则

21、执行(3)。(3) 如果ik¹k，则交换第k行与第ik行对应元素位置，。(4) 消元，对i = k, , n，计算lik = aik/ akk，对j = k + 1, , n+1，计算：aij = aijlikakj步骤 2：回代过程：(1) 若ann = 0，则矩阵奇异，程序结束；否则执行(2)。(2) xn = an,n+1 / ann，对i = n 1, , 2, 1，计算。五、实验结果与分析(一) 实验源程序function x=liezhuyuan(A,b) n=length(b);x=zeros(n,1);c=zeros(1,n);t=0;for i=1:n-1max=a

22、bs (A(i,i);m=i;for j=i+1:nif max<abs(A(j,i);m=j;endendif m=ifor k=1:nc(k)=A(i,k);A(i,k)=A(m,k);A(m,k)=c(k);endt=b(i);b(i)=b(m);b(m)=t;endfor k=i+1:nfor j=i+1:nA(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i);endb(k)=b(k)-b(i)*A(k,i)/A(i,i);A(k,i)=0;endendx (n) =b(n)/A(n,n);for i=n-1:-1:1sum=0;for j=i+1:nsum =s

23、um+A(i,j)*x(j);endx(i)=(b(i)-sum)/A(i,i);end(二) 实验结果>> A=0.101 2.304 3.555;-1.347 3.712 4.623;-2.835 1.072 5.643;b=1.183 2.137 3.035;x=liezhuyuan(A,b)x = -0.3982 0.0138 0.3351(二) 实验结果分析高斯消去法是通过选择选主元素，以达到更精确的数值，比原先没有用选主元素的方法得到的数值更精确。实验四解线性方程组的迭代法一、实验目的和意义1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较；2、运用所学的

24、迭代法算法，解决各类线性方程组，编出算法程序；3、体会上机计算时，终止步骤或k > (给予的迭代次数)，对迭代法敛散性的意义；4、体会初始解x(0)，松弛因子的选取，对计算结果的影响。二、实验内容对实验三所列目的和意义的线性方程组，试分别选用Jacobi 迭代法，Gauss-Seidol迭代法和SOR方法计算其解。三、实验设备和实验环境安装有C、C+或MATLAB的计算机。四、实验要求：1、 对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与改进平方根法；追赶法求解（选择其一）；2、 应用结构程序设计编出通用程序；3、 比较计算结果，分析数值解误差的原因；4

25、、 尽可能利用相应模块输出系数矩阵的三角分解式。实验五非线性方程求根一、实验目的和意义1、通过实验进一步了解方程求根的算法；2、认识选择计算格式的重要性；3、掌握迭代算法和精度控制；4、明确迭代收敛性与初值选取的关系。二、实验内容设方程f(x) = x3 3x 1 = 0有三个实根。现采用下面六种不同计算格式，求f(x) = 0的根或：1、2、3、4、5、6、三、实验设备和实验环境安装有C、C+或MATLAB的计算机。四、实验要求：1、编制一个程序进行运算，最后打印出每种迭代格式的敛散情况；2、用事后误差估计来控制迭代次数，并且打印出迭代的次数；3、初始值的选取对迭代收敛有何影响；4、分析迭代

26、收敛和发散的原因。实验六函数插值方法一、实验目的和意义1、 学会常用的插值方法，求函数的近似表达式，以解决其它实际问题；2、 明确插值多项式和分段插值多项式各自的优缺点；3、 熟悉插值方法的程序编制；4、 绘出插值函数的曲线，观察其光滑性。二、实验内容对于给定的一元函数y = f(x)的n+1个节点值yj = f(xj)，j = 0, 1, , n。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。数据如下：(1) 求五次Lagrange多项式L5(x)和分段三次插值多项式，计算f(0.596)，f(0.99)的值。(提示：结果为f(0.596) » 0.6

27、25732，f(0.99) » 1.05423)。xj0.40.550.650.800.951.05yj0.410750.578150.696750.901.001.25382(2) 试构造Lagrange多项式L6(x)，计算f(1.8)的值。结果f(1.8) » 0.164762，f(6.15) » 0.001266。xj1234567yj0.3680.1350.0500.0180.0070.0020.001三、实验设备和实验环境安装有C、C+或MATLAB的计算机。四、实验要求：1、利用Lagrange插值公式编写出插值多项式程序设；2、给出插值多项式或分段

28、三次插值多项式的表达式；3、根据节点选取原则，对问题(2)用三点插值或二点插值，观察其结果如何；4、对此插值问题如果用Newton插值多项式，其结果如何？Newton插值多项式如下：其中，。提高：多项式插值的振荡现象一、实验目的和意义在一个固定的区间上用插值逼近一个函数，显然Lagrange插值中使用的节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，Ln(x)是否也更加靠近被逼近的函数。Runge给出的一个例子是极著名并富有启发性的。二、实验内容设区间-1, 1上函数f(x)=1/(1+25x2)。考虑区间-1,1的一个等距划分，分点为xi= -1 + 2i/n，i=0，

29、1，2，n，则拉格朗日插值多项式为：。其中，li(x)，i=0，1，2，n是n次Lagrange插值基函数。三、实验设备和实验环境安装有C、C+或MATLAB的计算机。四、实验源程序function czzd% 数值实验：“多项式插值的震荡现象”% 输入：函数式选择，插值结点数% 输出：拟合函数及原函数的图形promps = '请选择实验函数，若选f(x),请输入f,若选h(x),请输入h,若选g(x),请输入g:'titles = 'czzd'result = inputdlg(promps,'czzd',1,'f');Nb_f

30、 = char(result);if(Nb_f = 'f' & Nb_f = 'h' & Nb_f = 'g')errordlg('实验函数选择错误！');return;endresult = inputdlg('请输入插值结点数N:','charpt_2',1,'10');Nd = str2num(char(result);if(Nd <1)errordlg('结点输入错误！');return;endswitch Nb_f case '

31、f' f=inline('1./(1+25*x.2)'); a = -1;b = 1; case 'h' f=inline('x./(1+x.4)'); a = -5; b = 5; case 'g' f=inline('atan(x)'); a = -5; b= 5;end x0 = linspace(a, b, Nd+1); y0 = feval(f, x0); x = a:0.1:b; y = Lagrange(x0, y0, x); fplot(f, a b, 'co'); hold

32、 on; plot(x, y, 'b-'); xlabel('x'); ylabel('y = f(x) o and y = Ln(x)-');%-function y=Lagrange(x0, y0, x);% Lagrange插值n= length(x0); m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if(j = k) p = p*(z - x0(j)/(x0(k) - x0(j); end end s = s + p*y0(k); end y(i) = s

33、;end五、结果分析(1) 增大分点n=2，3，时，拉格朗日插值函数曲线如图所示。n=3 n=6n=7 n=8n=9n=10从图中可以看出，随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f(x)，但是却在两端x= -1和x=1处出现了很大的振荡现象。并且，仔细分析图形，可以看出，当n为奇数时，虽然有振荡，但振荡的幅度不算太大，n为偶数时，其振荡幅度变得很大。通过思考分析，认为可能的原因是f(x)本身是偶函数，如果n为奇数，那么Lagrange插值函数Ln(x)的最高次项xn-1是偶次幂，比较符合f(x)本身是偶函数的性质；如果n为偶数，那么Lagrange插值函数Ln(x)的

34、最高次项xn-1是奇次幂，与f(x)本身是偶函数的性质相反，因此振荡可能更剧烈。(2) 将原来的f(x)换为其他函数如h(x)、g(x)，结果如图所示。其中h(x)，g(x)均定义在-5, 5区间上，h(x)=x/(1+x4)，g(x)=arctanx。h(x), n=7 h(x), n=8h(x), n=9 h(x), n=10g(x), n=8 g(x), n=9g(x), n=12 g(x), n=13分析两个函数的插值图形，可以看出：随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f(x)，但是却在两端x= -5和x=5处出现了很大的振荡现象。并且，仔细分析图形，可以看

35、出，当n为偶数时，虽然有振荡，但振荡的幅度不算太大，n为奇数时，其振荡幅度变得很大。原因和上面f(x)的插值类似，h(x)、g(x)本身是奇函数，如果n为偶数，那么Lagrange插值函数Ln(x)的最高次项xn-1是奇次幂，比较符合h(x)、g(x)本身是奇函数的性质；如果n为奇数，那么Lagrange插值函数Ln(x)的最高次项xn-1是偶次幂，与h(x)、g(x)本身是奇函数的性质相反，因此振荡可能更剧烈。扩展：样条插值的收敛一、实验目的和意义多项式插值是不收敛的，即插值的节点多，效果不一定好。样条函数插值能较好地解决这个问题，通过本实验可以验证这一理论结果。二、实验内容通过对前面实验中

36、的三个函数进行三次样条插值，与Lagrange插值进行对比，验证样条插值的收敛性。三、实验设备和实验环境安装有C、C+或MATLAB的计算机。四、实验源程序function ytcz% 数值实验： “样条插值的收敛”% 输入：函数式选择，插值结点数% 输出：拟合函数及原函数的图形promps = '请选择实验函数，若选f(x),请输入f,若选h(x),请输入h,若选g(x),请输入g:'titles = 'ytcz'result = inputdlg(promps,'ytcz',1,'f');Nb_f = char(result)

37、;if(Nb_f = 'f' & Nb_f = 'h' & Nb_f = 'g')errordlg('实验函数选择错误！');return;endresult = inputdlg('请输入插值结点数N:','charpt_2',1,'10');Nd = str2num(char(result);if(Nd <1)errordlg('结点输入错误！');return;endswitch Nb_f case 'f' f=inline

38、('1./(1+25*x.2)'); a = -1;b = 1; case 'h' f=inline('x./(1+x.4)'); a = -5; b = 5; case 'g' f=inline('atan(x)'); a = -5; b= 5;endx0 = linspace(a, b, Nd+1); y0 = feval(f, x0);x = a:0.1:b;cs = spline(x0, y0); y = ppval(cs, x);plot(x0, y0, 'o'); hold on; pl

39、ot(x, y, 'k-');xlabel('x'); ylabel('y = f(x) o and y = Spline(x)-');五、结果分析(1) 实验结果如图所示。f (x), n=5 f(x), n=10f (x), n=20 f(x), n=30h(x), n=5 h(x), n=10h(x), n=20 h(x), n=30g(x), n=5 g(x), n=10g(x), n=20 g(x), n=30从图中可以看出，三次样条插值随着插值结点的增加，由于其采用了分段三次多项式拟合的方法，因此并没有出现振荡现象。(2) 样条插值思

40、想最早产生于工业部门，如表，某汽车制造商用三次样条插值设计车门曲线，其中一段数据如表所示。xk012345678910yk0.00.791.532.192.713.033.272.893.063.193.29yk0.80.2Matlab源程序如下：function CarDoorDesignx0=0:10;y0=0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29;x=0:0.1:10;pp=csape(x0,y0,'complete',0.8 0.2);y = ppval(pp, x);plot(x0, y0, '

41、o'); hold on; plot(x, y, 'k-');xlabel('x'); ylabel('y = f(x) o and y = Spline(x)-');设计的车门曲线如图所示。实验七数值积分与数值微分一、实验目的和意义1、 深刻认识数值积分法的意义；2、 明确数值积分精度与步长的关系；3、 根据定积分的计算方法，可以考虑二重积分的计算问题。二、实验内容选用复合梯形公式，复合Simpson公式，Romberg算法，计算：(1) (2) (3) (4) 三、实验设备和实验环境安装有C、C+或MATLAB的计算机。四、实验要求：

42、1、 编制数值积分算法的程序；2、 分别用两种算法计算同一个积分，并比较其结果；3、分别取不同步长h = (b-a)/n，试比较计算结果(如n = 10, 20等)；4、 给定精度要求，试用变步长算法，确定最佳步长。实验八常微分方程初值问题数值解法一、实验目的和意义1、 熟悉各种初值问题的算法，编出算法程序；2、 明确各种算法的精度与所选步长有密切关系；3、 通过计算更加了解各种算法的优越性。二、实验内容科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler法，Rung-Kutta方法求其数值解，诸如以下问题：(1) 分别取h = 0.1, 0.2, 0.4时的数值解；

43、初值问题的精确解。(2) 取步长h = 0.1，用四阶标准R-K方法求值。(3) 用改进Euler法或四阶标准R-K方法求解取步长h = 0.01，计算y(0.05)，y(0.10)，y(0.15)数值解，参考结果y1(0.15) » -0.9880787，y2(0.15) » 0.1493359，y3(0.15) » 0.8613125。(4) 利用四阶标准R-K方法求二阶方程初值问题的数值解： (I) (II) (III) (IV) 三、实验设备和实验环境安装有C、C+或MATLAB的计算机。四、实验要求：1、根据初值问题数值算法，分别选择两个初值问题编程计算

44、；2、试分别取不同步长，考察某节点xj处数值解的误差变化情况；3、试用不同算法求解某初值问题，结果有何异常；4、分析各个算法的优缺点。实验九函数逼近与曲线拟合一、实验目的和意义1、掌握曲线拟合的最小二乘法；2、最小二乘法亦可用于解超定线代数方程组；3、探索拟合函数的选择与拟合精度间的关系。二、实验内容从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。t(分)0510152025303540455055y(´10-4)01.272.162

45、.863.443.874.154.374.514.584.024.64三、实验设备和实验环境安装有C、C+或MATLAB的计算机。四、实验要求：1、用最小二乘法进行曲线拟合；2、近似解析表达式为；3、打印出拟合函数，并打印出与的误差，j = 1, 2, , 12；4、另外选取一个近似表达式，尝试拟合效果的比较；5、绘制出曲线拟合图。附实验一：多项式最小二乘拟合一、实验目的编制以函数xkk=0,n；为基的多项式最小二乘拟合程序。二、实验内容对表中的数据作三次多项式最小二乘拟合。xi-1.0-0.50.00.51.01.52.0yi-4.447-0.4520.5510.048-0.4470.549

46、4.552取权函数wi1，求拟合曲线中的参数、平方误差2，并作离散数据xi, yi的拟合函数的图形。三、实验设备和实验环境安装有C、C+或MATLAB的计算机。四、实验源程序function charpt3% 数值实验：“zxecnh”% 输出：原函数及求得的相应插值多项式的函数的图像以及参数alph和误差rx0 = -1:0.5:2;y0 = -4.447 -0.452 0.551 0.048 -0.447 0.549 4.552;n = 3; % n为拟合阶次alph = polyfit(x0, y0, n);y = polyval(alph, x0);r = (y0 -y)*(y0 -y

47、)' % 平方误差x = -1:0.01:2;y = polyval(alph, x);plot(x, y, 'k-');xlabel('x'); ylabel('y0 * and polyfit.y-');hold onplot(x0, y0, '*')grid on;disp('平方误差:', num2str(r)disp('参数alph:', num2str(alph)五、实验结果分析平方误差：2.1762e-005参数alph：1.9991 -2.9977 -3.9683e-005

48、0.54912附实验二：正交化多项式的最小二乘拟合一、实验目的编制正交化多项式最小二乘拟合程序。二、实验内容用正交化多项式最小二乘拟合求解附实验一中的三次多项式最小二乘拟合问题，作拟合曲线的图形，计算平方误差。三、实验设备和实验环境安装有C、C+或MATLAB的计算机。四、实验源程序function charpt3% 数值实验:“实验zjhdxszxecnh”% 子函数调用:dlsa% 输出:原函数及求得的相应插值多项式的函数的图像以及参数alph和误差rx0 = -1:0.5:2;y0 = -4.447 -0.452 0.551 0.048 -0.447 0.549 4.552; resul

49、t = inputdlg('请输入权矢量w:','charpt_3', 1, '1 1 1 1 1 1 1'); w = str2num(char(result); n = 3; % n为拟合阶次 a, b, c, alph, r = dlsa(x0, y0, w, n);disp('平方误差:', num2str(r)disp('参数alph:', num2str(alph)%-function a, b, c, alph, r = dlsa(x, y, w, n)% 功能：用正交化方法对离散数据作多项式最小二乘

50、拟合m = length(x) - 1;if(n < 1| n >= m)errordlg('错误:n<1或者n>=m!');return;end% 求三项递推公式的参数a, b, 拟合多项式s(x)的系数c,其中d(k) = (y, sk);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'% 求拟合多项式s(x)的降幂系数alphalph=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: