《数值分析》课程考试

1、淮 阴 工 学 院数值分析考试基于Matlab的方法综合应用报告班 级： 金融1121 姓 名： 蒋倩 学 号： 1124104119 成 绩： 数 理 学 院 2014年6月7日数值分析课程考试2014.6.3一 课程主要教学内容综述（开卷部分）（70分）综述要求：1. 综述内容以本课程实际讲授内容和所使用的教材为基准；2. 按主题论述，问题提出的背景、动因，解决问题的理论基础和方法；3. 以数值计算的实例说明方法的应用，并进一步对方法本身和应用结果进行简要的分析和评价；4. 各数值计算的Matlab程序按顺序编号，附在正文的后面。*正文用五号字，附件编号，用小五号字，至少五页，可以讨论，要

2、用自己的语言*参考文献：书（白峰杉，数值计算引论（第二版），北京：高等教育出版社，2010）、电子讲义（章节名，电子讲义，淮阴工学院，2014）*自己再找一到两个参考文献*格式：正文、参考文献、附录正 文一、线性方程组求解的数值方法1、高斯消去法：1.1问题提出的背景、动因1：高斯消去法是求解线性代数方程组的直接方法，假设计算中没有舍入误差，经过有限次算数运算能够给出问题的精确解。1.2解决问题的理论基础和方法1,2：高斯消去法的基本思想就是将一般的线性方程组化成与之等价的上三角或下三角形式来求解。高斯消去法的消元过程，从代数运算的角度看就是用一个下三角矩阵左乘方程组的系数矩阵，且乘积的结果为

3、上三角矩阵，即LU 分解的MATLAB实现： 或 1.3例：Ax=b,已知A = 1 -2 0.33 0;2 -4.01 -1 5;-1 3 1 4; -3 6 0 2，b = -1 1 0 3'，用高斯消去法求解x（程序见附录一）结果：X = -3.70 -1.47 -0.72 0.36L = 1.00 0 0 0 0.33 1.00 0 0 -0.67 -0.01 1.00 0 -0.33 0 -0.33 1.00U = -3.00 6.00 0 2.00 0 1.00 1.00 3.33 0 0 -0.99 6.37 0 0 0 2.79P = 0 0 0 1.00 0 0 1.

4、00 0 0 1.00 0 0 1.00 0 0 0x = -3.70 -1.47 -0.72 0.361.4分析和评价：用高斯消去法(LU分解)和用Matlab求解的结果相同，所以这种方法较为适用。2、迭代法2.1问题提出的背景、动因3：通常逆矩阵不易求得，特别是对于大型的线性方程组，需要用迭代法求解。2.2解决问题的理论基础和方法3：用迭代法求解线性方程组，首先要把线性方程组写成等价的形式(1)由迭代格式（1）确定如下的迭代算法：Jacobian 迭代法:Gauss-Seidel 迭代法:2.3例：生成一个矩阵A和一个向量b，用Gauss迭代法和Jacobi迭代法求解方程组。（程序见附录二

5、）结果：X0 = 0.86 0.19 -0.20 0.32 -0.13IterN0 = 16.00X1 = 0.86 0.19 -0.20 0.32 -0.13IterN1 = 40.002.4分析和评价：Jacobian迭代法与Gass-seidel迭代,迭代次数比较：Jacobian迭代法在要求的精度下求得近似解需要的次数为16次，Gass-seidel迭代法在要求的精度下求得近似解需要的次数为40次。Gass-seidel迭代法迭代的次数更加的少，即Gass-seidel迭代法在迭代次数上比Jacobian迭代法更加的优越。二、范数1、问题提出的背景、动因1：范数是用于度量“量”大小的概

6、念，是进行算法分析的基本工具，在数值计算中对方法的评价，几乎总是离不开范数。2、解决问题的理论基础和方法9：向量的范数：p-范数 矩阵（算子）的范数3、例：生成一个5阶矩阵，并计算它的1-范数，2-范数和无穷范数。（程序见附录三）结果：M的1范数为： 5.65M的2范数为： 3.75M的无穷范数为： 6.144、分析和评价：用norm函数，实现对1、2、无穷范数的求解，和用公式计算法相比方法简单，易于实现。三、插值1、Lagrange插值1.1问题提出的背景、动因5：在科学研究和计算中，往往回遇到复杂函数的分析与计算，有时用简单的函数来代替，可能会去掉不必要的麻烦而使问题比较容易地得到解决。1

7、.2解决问题的理论基础和方法5：构造插值多项式的基函数：， 因为拉格朗日插值多项式的基函数有如下的性质：，所以拉格朗日插值多项式， 满足插值的条件。1.3例：给定函数：，. （1） 借助Matlab 求该函数的Lagrange 插值基函数以及差值多项式的表达式。（程序见附录四）L = (4398046511104*(pi - x)*(pi - 2*x)*(pi - 4*x)*(3*pi - 4*x)/1285229138915719, (70368744177664*x*(pi - x)*(pi - 2*x)*(3*pi - 4*x)/1285229138915719, -(281474976

8、710656*x*(pi - x)*(pi - 4*x)*(3*pi - 4*x)/6854555407550501, (70368744177664*x*(pi - x)*(pi - 2*x)*(pi - 4*x)/1285229138915719, -(4398046511104*x*(pi - 2*x)*(pi - 4*x)*(3*pi - 4*x)/1285229138915719Ln = (35184372088832*2(1/2)*x*exp(-(9*pi2)/80)*(9*pi2)/16 + 1)(1/2)*(pi - x)*(pi - 2*x)*(pi - 4*x)/12852

9、29138915719 - (4398046511104*x*exp(-pi2/5)*(pi2 + 1)(1/2)*(pi - 2*x)*(pi - 4*x)*(3*pi - 4*x)/1285229138915719 - (4398046511104*(pi - x)*(pi - 2*x)*(pi - 4*x)*(3*pi - 4*x)/1285229138915719 - (35184372088832*2(1/2)*x*exp(-pi2/80)*(pi2/16 + 1)(1/2)*(pi - x)*(pi - 2*x)*(3*pi - 4*x)/12852291389157191.4分析

10、和评价：对于较简单的函数而言，用公式即可求出Lagrange 插值基函数以及差值多项式，而对于较复杂的函数用Matlab计算易于实现，且结果准确。2、三次样条插值2.1问题提出的背景、动因1：随着x的值越来越大，lagrang模拟的情况会越来越差，原因就是lagrang插值法出现了“龙格”现象。因此，当x值较大时，可以选用三次样条插值。2.2解决问题的理论基础和方法5：通常在插值区间的端点处附加2个条件：第一类边界条件：固定端点的斜率：固定边界条件：， 第二类边界条件：给定端点的的二阶导数：自由边界条件：，. 第三类边界条件：周期性条件：2.3例：取一个复杂方程，设置一阶边界条件、二阶边界条件

11、以及混合边界条件。（程序见附录五）结果：一阶边界条件：DY = 0.5000 0.0000D2Y = -12.7199 - 0.1857i -0.0067 + 0.0016iDy = -0.8660 0.0000 + 0.0006iD2y = -2.0000 -0.0000 - 0.0040i二阶边界条件：DY = -1.4071 - 0.0241i 0.2604 - 0.0002iD2Y = 2.0000 2.0000 + 0.0000iDy = -0.8660 0.0000 + 0.0006iD2y = -2.0000 -0.0000 - 0.0040i混合边界条件：DY = 2.0000

12、 0.1301 - 0.0002iD2Y = -24.2980 - 0.1857i 1.0000 + 0.0000iDy = -0.8660 0.0000 + 0.0006iD2y = -2.0000 -0.0000 - 0.0040i2.4分析和评价：对于同一个函数，根据两端端点的倒数设置不同的边界条件导致拟合的效果不同。四、最小二乘拟合1、问题提出的背景、动因10：在生产实践中，人们常常需要根据已知观测数据确定不同量之间的关系，为进一步判断、预测等提供理论依据，数据拟合的最小二乘法正是解决这个问题的重要方法。2、解决问题的理论基础和方法6：当拟合多项式用一组多项式的线性组合时 到距离的平方

13、是组合系数的元二次函数： 所以最小二乘逼近多项式必须满足如下必要条件： 即满足法方程组: 3、例：XI = 0 0.8 1.6 2.8 3.2 4.8.5.2 6.4 7.2 8.0，YI = 1.00 1.49 2.23 4.32. 5.55 7.97 7.80 7.20. 7.08 7.84，用最小二乘法拟合该组数据。（程序见附录六）结果：n=5时n=9时4、分析和评价：当有n个数据进行拟合时，最高拟合次数不得超过n-1次，而且拟合次数越大，拟合效果则越高，精确度也就越高。五、数值积分1、两点公式（梯形公式）和三点公式（Simpson公式）1.1问题提出的背景、动因1：随着符号计算方法及相

14、应数学软件工具的发展，计算积分又有了一种有力途径，但众多的实际问题中，积分计算更多的还要依靠数值方法。经典的数值积分方法都可以看成是从几何直观出发构造的。1.2解决问题的理论基础和方法7：两点公式（梯形公式）： 2. 三点公式（Simpson公式）： 1.3例：梯形求积公式和Simpson求积公式在积分区间上分别用直线和抛物线 代替曲线。（程序见附录七：程序1：Quad_Example_01.m、程序2：Quad_Example_02.m）结果：图1图21.4分析和评价：由图1和图2还表明，在整个积分区间上用直线或用抛物线代替原曲线做积分，尽管计算简单，但精度较差，结果不好。2、复化求积公式2

15、.1问题提出的背景、动因7：在整个积分区间上用直线或用抛物线代替原曲线做积分，尽管计算简单，但精度较差，结果不好，因此，一个自然的想法是：在尽可能保持公式计算简单、易于计算机实现的优点的前提下，寻找提高误差精度的途径和方法，即复化求积公式。2.2解决问题的理论基础和方法7：复化的梯形求积公式 复化的Simpson求积公式 2.3例：X0=0 0.20 1.00 2.10 3.50 5.00 6.80 7.50 9.00 11.2 12.0，Y0=-1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86，试利用复化的梯形求积法和复化的Simp

16、son求积法求该组数据所在曲线与基准线（轴）在范围内所围成图形面积.画出数据散点图和图形的示意图。（程序见附录八）结果：梯形求积公式给出的面积： 11.6090simpson求积公式给出的面积为： 11.91282.4分析和评价：两种方法都易于计算机实现，Simpson求积法所产生的边缘曲线具有一定的光滑性（一阶导数连续），这是复化梯形求积法所不具有的优点。3、Monte Carlo法3.1问题提出的背景、动因1：当空间维数大于1时，按前面介绍的数值方法的思想是很难处理的，此时可以运用Monte Carlo方法。3.2解决问题的理论基础和方法7：函数在区间上的平均值为当，那么有(31)由此，(

17、32)利用在上服从均匀分布的随机数生成函数unifrnd(a,b,1,n)生成，代入（32）3.3例：X0=0 0.20 1.00 2.10 3.50 5.00 6.80 7.50 9.00 11.2 12.0，Y0=1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86，利用Monte Carlo法求面积问题。（程序见附录九）结果：MC_ls = 11.4987MC_spline = 11.55263.4分析和评价：由于Monte Carlo 法是利用随机函数生成的随机数求函数的平均值，所以每次计算的结果会有所不同。六、非线性代数方程（组

18、）的数值解法1、二分法1.1问题提出的背景、动因10：对于无规律的非线性方程的求解也无精确解法，因此，研究非线性方程的数值解法成为必然。求根方法中最直观最简单的方法是二分法。1.2解决问题的理论基础和方法4：二分法的基本思想是通过计算隔根区间的中点，逐步将隔根区间缩小，从而可得方程的近似根数列。二分法的求根过程：用表示方程在区间上的精确解，对于给定的精度要求，取区间的中点，并按下式进行判断： 1.3例：用二分法确定函数的根。（程序见附录十：Find_Interval_of_Including_Root_01.m）图像：1.4分析和评价：二分法适用于一个方程的场合，收敛速度是线性的。用二分法求解

19、一元方程较为适用，易于实现。但求解二元方程就比较复杂了，所以二分法有局限性。2、 黄金分割法2.1问题提出的背景、动因10：黄金分割法和二分法都是求非线性方程的有效手段，它们本质相同，不同点在于黄金分割法取区间的黄金分割点，而二分法只是取中点而已。2.2解决问题的理论基础和方法4：在区间内取对称的两点： 使得求根（方法）程序如下： 2.3例：用黄金分割法求的根。（程序见附录十一：Gold_Rate_method_01.m）结果：X = 2.22182.4分析和评价：由于黄金分割法把区间分割为不等的两个区间，因此黄金分割法优于二分法，且其达到要求的精确度所需计算的次数小于二分法的计算次数。3、

20、牛顿迭代法3.1问题提出的背景、动因10：二分法在解一元方程应用中简单，特别是写程序要简单些，但是求根过程的步骤较多。为加快求根速度，可以使用牛顿迭代法。3.2解决问题的理论基础和方法4：牛顿迭代格式： 以及迭代法: 3.3例：.求方程 的实根。 （程序见附录十二：Newton_Iter_Fun_01.m、Newton_Iteration_method_03.m） 结果： 在-0.7附近的实根：在1.2附近的实根：3.4分析和评价：牛顿迭代法不仅可以解一元方程，还可以解多元可微分方程，而且求根速度很快。二分法和黄金分割法都是1阶收敛的。Newton 迭代法至少是2阶收敛的。参 考 文 献参考文

21、献1书（白峰杉，数值计算引论（第二版），北京：高等教育出版社，2010）2电子讲义（Lecture2_线性方程组，电子讲义，淮阴工学院，2014）3电子讲义（Lecture3_解线性方程组的迭代法，电子讲义，淮阴工学院，2014）4电子讲义（Lecture4_非线性代数方程（组）的数值解法，电子讲义，淮阴工学院，2014）5电子讲义（Lecture5_插值，电子讲义，淮阴工学院，2014）6电子讲义（Lecture6_最小二乘拟合与最佳逼近，电子讲义，淮阴工学院，2014）7电子讲义（Lecture7_微积分的数值方法，电子讲义，淮阴工学院，2014）8电子讲义（Lecture8_微分方程（组

22、）初值问题数值方法，电子讲义，淮阴工学院，2014）9电子讲义（内积与范数，电子讲义，淮阴工学院，2014）10网上资料附 录附录一：% LU Factorization Method%clear allclcformat bankA = 1 -2 0.33 0; 2 -4.01 -1 5;-1 3 1 4; -3 6 0 2b = -1 1 0 3'X = Abn = length(b);Rx = 1e-6;% Solve the Equations by LU Factorization Method% L*U = P*AL U P = lu(A)b = P*b;y(1) = b(

23、1);for i = 2:n Sum = 0; for j = 1:i -1 Sum = Sum + L(i, j)*y(j); end y(i) = b(i) - Sum;endx(n) = U(n, n)y(n);for i = n-1:-1:1 Sum = 0; for j = i + 1 : n Sum = Sum + U(i, j)*x(j); end x(i) = U(i, i)(y(i) - Sum);end x = x'附录二：clearclcn = 5;A = randn(n);for i = 1:n A(i,i) = sum(abs(A(i,:);endb = ra

24、ndn(n, 1);D = diag(diag(A);L = -tril(A, -1);U = -triu(A,1); % Gauss-Seidel Iteration MethodM0 = (D - L)U;f0 = (D - L)b;Rho0 = max(abs(eig(M0)R = 1e-08;switch sign(1 - Rho0) case -1 disp('The Gauss-Seidel Method is not applicable') case 0 disp('The Gauss-Seidel Method is not applicable

25、9;) otherwise x(:,1) = normrnd(0, 9, n, 1); k = 1; while k<=50*n x(:, k+1) = M0*x(:, k) + f0; if norm(x(:, k+1) - x(:, k)>=R k = k+1; else X0 = x(:, k+1) IterN0 = k break end endend % Jacobian Iteration MethodM1 = D(L+U);f1 = Db;Rho1 = max(abs(eig(M1)R = 1e-08;switch sign(1 - Rho1) case -1 dis

26、p('The Jacobian Method is not applicable') case 0 disp('The Jacobian Method is not applicable') otherwise x(:,1) = normrnd(0, 9, n, 1); k = 1; while k <= 50*n x(:, k+1) = M1*x(:, k) + f1; if norm(x(:, k+1) - x(:, k)>=R k = k+1; else X1 = x(:, k+1) IterN1 = k break end endend附录三

27、：clear allclcM=normrnd(0,1,5);Norm1=norm(M,1); Norm2=norm(M,2); Normi=norm(M,inf); disp(' M的1范数为：') disp(Norm1)disp(' M的2范数为：')disp(Norm2)disp(' M的无穷范数为：')disp(Normi)附录四：%The Basic Function of Lagrange Interpolationclear allclcsyms x t f Lx = sym(x, 'real');t = sym(0

28、: pi/4 : pi);y = sqrt(1+t.2).*sin(t-pi/2).*exp(-t.2/5);A = ones(length(t), 1); B = factor(prod(x - t);N = 1:length(t);for k = N for i = N if i = k A(k) = A(k)*(t(k) - t(i); else A(k) = A(k); end end L(k) = B/(x - t(k)*A(k);endL = simple(L)Ln = sum(L.*y)附录五：clear clcclf% Use for Plotting Curve of the

29、 interpolated functiont = sym('t');f = sqrt(1 - t2)*cos(t + pi/3)*exp(-t2);Df = simple(diff(f);D2f = simple(diff(f, 2);fun = matlabFunction(f);Dfun = matlabFunction(Df);D2fun = matlabFunction(D2f); % the Interpoaltion NodesXI = linspace(0, pi, 8);YI = fun(XI);DY = Dfun(XI(1), XI(end);D2Y = D

30、2fun(XI(1), XI(end); x = linspace(0, pi, 51);y = fun(x);set(gca,'FontSize', 24) plot(x, y, 'r-', 'linewidth', 3, 'markersize', 10)pause(3)hold on% PP = csape(XI, YI, 1, 1, 0.5, 0); %一阶边界条件% PP = csape(XI, YI,2.2, 1, 0.5); %二阶边界条件% PP = csape(XI, YI, 1, 2, 2, 1) %混合边界条

31、件 % linewidth = 3, The interval is a, bpause(3)fnplt(PP, 3, -0.1, pi + 0.1)xlim(-0.1, pi + 0.1) Fun = fnval(PP, x);F = ppval(PP, x);% plot(x, F, 'linewidth', 3, 'markersize', 10) plot(XI, YI, 'ko', 'linewidth', 3, 'markersize', 10)% legend('The Curve of Th

32、e Function', 'Cubic Spline', 'Data')title('The Cubic Spline Interpolation With End Conditions')% The derivatives at ends of the spline functionDY = fnval(fnder(PP), x(1), x(end) D2Y = fnval(fnder(PP, 2), x(1), x(end) % The derivatives at ends of the interpolantDf = diff(f

33、, t);D2f = diff(f, t, 2);Dy = subs(Df, t, x(1), x(end)D2y = subs(D2f, t, x(1), x(end)% The Coeficients of the interpolation polynomialspause(3)I = 0;for k = 1 : length(XI) - 1 SP = poly2sym(PP.coefs(k, :); P = matlabFunction(SP); X=linspace(XI(k),XI(k+1),21); plot(X,P(X-XI(k),'k','linewi

34、dth',3)endSI = sum(abs(I)附录六：clearclcXI = 0 0.8 1.6 2.8 3.2 4.8.5.2 6.4 7.2 8.0;YI = 1.00 1.49 2.23 4.32. 5.55 7.97 7.80 7.20 7.08 7.84; n = 5;% n = 9;P, S = polyfit(XI, YI, n)t = linspace(XI(1), XI(end), 41);% Use the fitting polynomial P to evaluationPP = polyval(P, t); plot(t, PP, XI, YI, 

35、9;ro', 'LineWidth', 3, 'markersize', 8)set(gca, 'FontSize', 24)legend('The Fitting Curve', 'The Data', 4)title('Curve Fitting by Least Square Approximation')附录七：程序1：% example of Quadrature Method quadclearclc Fun = (x) x.*cos(2*x) + 2;x = linspace(

36、-pi/4, pi/3, 61);f = Fun(x);y0 = 0: f(1)/20 : f(1);ye = 0: f(end)/20 : f(end);ylin = (x(end) - x(1)(f(end) - f(1)*(x - x(1) + f(1);set(gca,'FontSize',24)plot(x, f, x, ylin, 'r', x(1) + 0*y0, y0, 'k', x(end) + 0*ye, ye, 'k', 'LineWidth', 3, 'markersize'

37、, 8)title('Numerical Integration by Trapezoid Quadrature', 'fontsize', 24)% legend('The Cuver of y = f(x)', 'The Secant Line')text(x(1), f(1) + 0.1, '(a, f(a)', 'fontsize', 24)text(x(end) + 0.05, f(end) + 0.05, '(b, f(b)', 'fontsize', 2

38、4)axis equalxlim(x(1) - 0.5, x(end) + 0.5);ylim(0, max(f) + 0.5);xlabel('X')ylabel('f(x)')程序2：% example of Quadrature Method quadclearclc Fun = (x) x.*cos(2*x) + 2;x = linspace(-pi/4, 3*pi/9, 61);f = Fun(x);y0 = 0: f(1)/20 : f(1);ye = 0: f(end)/20 : f(end);n = 2;h = (x(end) - x(1)/n;

39、XI = x(1), x(1) + h, x(end);YI = Fun(x(1), Fun(x(1) + h), Fun(x(end);yparabola = Lagrange_Fun_01(x, XI, YI);set(gca,'FontSize',24)plot(x, f, x, yparabola, 'r', x(1) + 0*y0, y0, 'k', x(end) + 0*ye, ye, 'k', 'LineWidth', 3, 'markersize', 8)title('Num

40、erical Integration By Simpson Quadrature', 'fontsize', 24)% legend('The Cuver of y = f(x)', 'The Parabola')text(x(1) - 0.4, f(1) + 0.1, '(a, f(a)', 'fontsize', 24)% text(x(1) + h - 0.5, Fun(x(1) + h) + 0.3, '(a + h, f(a + h)', 'fontsize', 2

41、4)text(x(end) + 0.05, f(end) + 0.1, '(b, f(b)', 'fontsize', 24)axis equalxlim(x(1) - 0.5, x(end) + 0.5);ylim(0, max(f) + 0.5);xlabel('X')ylabel('f(x)') 附录八：程序1：clearclcX0=0 0.20 1.00 2.10 3.50 5.00 6.80 7.50 9.00 11.2 12.0;Y0=-1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0

42、.39 0.77 0.86;T=0; for i=2:length(X0) T=T+(Y0(i)+Y0(i-1)*(X0(i)-X0(i-1)/2; enddisp('梯形求积公式给出的面积：')disp(-T)figure(1)set(gca,'fontsize',14)patch(X0(1) X0 X0(end),0 Y0 0,'r')hold onplot(X0,Y0,'k.',X0,Y0,'b','markersize',20,'linewidth',2.5)plot(X0,

43、0.*X0,'k','linewidth',3)ylim(-2,0.2)程序2：clearclcX0=0 0.20 1.00 2.10 3.50 5.00 6.80 7.50 9.00 11.2 12.0;Y0=-1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86;T = 0;for k = 1: 2: length(X0) - 2 T = T + 6(X0(k+2) - X0(k)*(Y0(k) + 4*Y0(k+1) + Y0(k+2); X(k, :) = linspace(X0(k), X0(k+

44、2), 20); Y(k, :) = Lagrange_Fun_01(X(k,:), X0(k : k+2), Y0(k : k+2);enddisp('simpson求积公式给出的面积：')disp(-T)figure(2)set(gca,'FontSize', 24)hold onfor k = 1: 2 :length(X0) - 2 patch(X0(k), X(k,:), X0(k +2), 0, Y(k, :), 0,'r') plot(X(k, :), Y(k, :), 'LineWidth', 3)endplot(

45、X0, Y0, 'o', X0, 0*X0, 'k', 'LineWidth', 3, 'markersize', 8)附录九：clcclearX0=0 0.20 1.00 2.10 3.50 5.00 6.80 7.50 9.00 11.2 12.0;Y0=1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86;Y_poly=polyfit(X0,Y0,8); X_data=0:0.01:pi; Y_data=polyval(Y_poly,X_data);Di = (Y0(2

46、) - Y0(1)/(X0(2) - X0(1);De = (Y0(end) - Y0(end - 1)/(X0(end) - X0(end - 1);PP = csape(X0, Y0, 'complete', Di, De); x=0:0.01:pi; y_spline = ppval(PP, x);N = 5000;L = 100;X = zeros(L, N);F_spline = (x)ppval(PP, x);for k = 1 : L X(k, :) = unifrnd(X0(1), X0(end), 1, N); MC_ls(k) = (X0(end) - X0

47、(1)*mean(polyval(Y_poly, X(k, :); MC_spline(k) = (X0(end) - X0(1)*mean(F_spline(X(k, :);endMC_ls = mean(MC_ls)MC_spline = mean(MC_spline)附录十：% Solving algebraic equations% Find an interval of including unque root%clear allclc% The equation: (x.5 - 2*x.3 +1).*sin(x) - 3*x.2 + 1.5 = 0;x = -2.5:0.05:3.

48、1;Fun_1 = (x) (x.5 - 2*x.3 +1).*sin(x) - 3*x.2 + 1.5;fun_1 = Fun_1(x);Fun_2 = (x) x.3 - x.2 -0.8*x +0.75;fun_2 = Fun_2(x);figure(1)plot(x, fun_1, x, 0*x, 'LineWidth', 3)title('Figure of f(x) = (x.5 - 2*x.3 +1).*sin(x) - 3*x.2 + 1.5', 'fontsize', 24)pause(2)% Find the interval

49、 of including unique root of the algebraic equationa = -3;% a = -1;% a = 1;% a = 2;Fa = Fun_1(a);Delta = 0.05;M = 30;n = 0;while n < M s = a + Delta; Fb = Fun_1(s); if Fa*Fb < 0 b = s; break else a = s; n = n + 1; endendt = linspace(a, b, 101);F = Fun_1(t);figure(2)plot(t, F, t, 0*t, 'LineWidth', 3)title('Figure of f(x) = (x.5 - 2*x.3 +1).*sin(x) - 3*x.2 + 1.5', 'fontsize', 24) 附录十二：function y = Newton_Iter_Fun_01(t)t2 = t.2;t3 = t2.2;y = t2.*(-3.0./2.0)-t3.*2.0+t.*t3+2.0;% Solving algebraic equations% Newton