数值计算课程7090307龙凯翔

1、东北大学秦皇岛分校数值计算课程设计报告系 别信息与计算科学专 业学 号7090307姓 名龙凯翔指导教师张建波 姜玉山成 绩教师评语：指导教师签字： 2011年7月13日1 绪 论数值分析是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。MATLAB是主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案。2 用M

2、ATLAB进行插值计算题目下列数据点的插值X01491625364964Y012345678可以得到平方根函数的近似，在区间0, 64上作图。(1) 用这9个点作8次多项式插值.(2) 用三次样条（第一边界条件）程序求.从得到的结果看在0, 64上，哪个插值更精确；在区间0, 1上，两种插值哪个更精确？2.1 程序文件建立新的m文件unit2.m：代码如下：clc; clear; x = 0, 1, 4, 9, 16, 25, 36, 49, 64; y = 0:8; xi = 0:64; p = polyfit(x, y, 8); yi = polyval(p, xi); subplot(2

3、, 2, 1); plot(x, y, 'x', xi, yi, 'k')xlabel('x轴'); ylabel('y轴'); legend('插值点', '8次差值多项式')yi = interp1(x, y, xi, 'cubic'); subplot(2, 2, 2); plot(x, y, 'x', xi, yi, 'k')xlabel('x轴'); ylabel('y轴'); legend('插值点&

4、#39;, '三次样条插值')xi = 0:0.01:1; p = polyfit(x, y, 8); yi = polyval(p, xi); subplot(2, 2, 3); plot(xi, yi, 'k')xlabel('x轴'); ylabel('y轴'); legend('8次差值多项式')yi = interp1(x, y, xi, 'cubic'); subplot(2, 2, 4); plot(xi, yi, 'k')xlabel('x轴'); y

5、label('y轴'); legend('三次样条插值') 图12.2 图样运行文件得到图12.3 分析比较两个函数的图像可知，在区间0, 64上三次样条插值函数要更精确一些，在区间0, 1上拉格朗日插值函数仍然不如三次样条插值函数更精确。3 多项式曲线拟合题目由实验给出数据表X0.00.10.20.30.50.81.0Y1.00.410.500.610.912.022.46试求3次、4次多项式的曲线拟合，再根据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。3.1程序文件建立新的m文件unit3.m：代码如下：clc; clear; x

6、 = 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; a = polyfit(x, y, 3); b = polyfit(x, y, 4); xi = 0.0:0.01:1.0; aa = polyval(a, xi); bb = polyval(b, xi); subplot(1, 2, 1); plot(x, y, 'x', xi, aa, 'k')xlabel('x轴'); ylabel('y轴'); legend

7、('插值点', '3次曲线拟合')subplot(1, 2, 2); plot(x, y, 'x', xi, bb, 'k')xlabel('x轴'); ylabel('y轴'); legend('插值点', '4次曲线拟合')poly2str(a, 'x')poly2str(b, 'x') 3.2 图形和结果运行程序结果如下：ans = - 6.6221 x3 + 12.8147 x2 - 4.6591 x + 0.92659ans

8、= 2.8853 x4 - 12.3348 x3 + 16.2747 x2 - 5.2987 x + 0.94272图像如下：图23.2 5次多项式曲线拟合下面进行的：代码如下：clc; clear; 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; a = polyfit(x, y, 5); xi = 0.0:0.01:1.0; aa = polyval(a, xi); plot(x, y, 'x', xi, aa, 'k')xlabel(&#

9、39;x轴'); ylabel('y轴'); legend('插值点', '5次曲线拟合')poly2str(a, 'x')运行程序结果如下：ans = - 79.3261 x5 + 195.4554 x4 - 172.7104 x3 + 69.0498 x2 - 11.0044 x + 0.99547图像如下：图3由图像可知，次数越高，拟合越准确。4 数值积分题目用不同数值方法计算积分.(1)取不同的步长h.分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小

10、的h，使得精度不能呢个再被改善？(2)用龙贝格求积计算完成问题(1).(3)用自适应辛普森积分，使其精度达到4.1 复合梯形公式代码建立新的m文件unit41.m：代码如下：clc, clear; x = 0.0000001:0.1:1; y = x.(1 / 2). * log(x); Fx = trapz(x, y)error = Fx + 4 / 9运行文件结果为：Fx = - 0.4123error = 0.0321取不同步长h:h = 0.01，Fx = - 0.4431，error = - 0.0014；h = 0.001，Fx = - 0.4444，error = - 5.487

11、5e - 005；h = 0.0001，Fx = - 0.4444，error = - 2.0338e - 006；h = 0.00001，Fx = - 0.4444，error = - 6.4496e - 008.可见，没有一个最小的h使精度不能再改变。4.2 复合辛普森公式代码建立新的m文件unit42.m：代码如下：Fx = quad(inline('x.(1 / 2). * log(x)'), 0.000001, 1，0.0001)error = Fx + 4 / 9运行文件结果为：Fx = - 0.4440error = 4.3273e - 004取不同的步长h：h

12、= 0.00001，Fx = - 0.4444，error = - 6.3537e - 005；h = 0.000001，Fx = - 0.4444，error = - 2.9938e - 006；h = 0.0000001，Fx = - 0.4444，error = - 3.0657e - 007；h = 0.0000000000001，Fx = - 0.4444，error = - 9.6548e - 009；h = 0.00000000000001，Fx = - 0.4444，error = - 9.6548e - 009。可见，存在一个最小的h使精度不能在改善，且h在0.0000000

13、000001与0.00000000000001之间。4.3龙贝格算法建立新的m文件unit43.m：代码如下：function q, n = unit43(f, a, b); M = 1; abs0 = 10; k = 0; T = zeros(1, 1); h = b - a; T(1, 1) = (h / 2) * (subs(f, a) + subs(f, b); while abs0>0.0001 k = k + 1; h = h / 2; p = 0; for i = 1:M x = a + h * (2 * i - 1); p = p + subs(f, x); end T(

14、k + 1, 1) = T(k, 1) / 2 + h * p; M = 2 * M; for j = 1:k T(k + 1, j + 1) = (4j) * T(k + 1, j) - T(k, j) / (4j - 1); end abs0 = abs(T(k + 1, j + 1) - T(k, j); endq = T(k + 1, k + 1); n = k; 在命令窗口输入：Fx, n = unit43('sqrt(x) * log(x)', 10( - 8), 1)输出结果为：Fx = - 0.4444n = 94.4 自适应辛普森算法Fx = quad(inl

15、ine('x.(1 / 2). * log(x)'), 0.000001, 1)输出结果为：Fx = - 0.4444。5 解线性方程题目线性方程组Ax = b的A及b为A = ，b = ，则解x = .用MATLAB内部函数求detA及A的所有特征值和cond(A) .若令A + A = ，求解(A + A)(x + x) = b, 输出向量x和.从理论结果和实际计算两方面分析线性方程组Ax = b解的相对误差及A的相对误差的关系.5.1 求detA及A的所有特征值和cond(A)输入程序：A = 10 7 8 7; 7 5 6 5; 8 6 10 9; 7 5 9 10de

16、t(A)V, D = eig(A)norm(A, 2)输出结果为：A = 10 7 8 7 7 5 6 5 8 6 10 9 7 5 9 10ans = 1V = 0.5016 - 0.3017 0.6149 0.5286 - 0.8304 0.0933 0.3963 0.3803 0.2086 0.7603 - 0.2716 0.5520 - 0.1237 - 0.5676 - 0.6254 0.5209D = 0.0102 0 0 0 0 0.8431 0 0 0 0 3.8581 0 0 0 0 30.2887ans = 30.2887则得到detA = 1A的特征值为0.0101 0.

17、8431 3.8581 30.2887 = 30.28876 解线性方程组的迭代法题目给出线性方程组其中系数矩阵为希尔伯矩阵：，i = 1, 2n.假设，若取n = 6, 8, 10，分别用雅克比迭代及SOR迭代(w = 1, 1.25, 1.5)求解.比较计算结果.6.1 用雅可比迭代法求解建立新的m文件unit61.m：代码如下：n = 6; H = hilb(n); b = H * ones(n, 1); e = 0.00001; for i = 1:n if H(i, i) = = 0 '对角元为零，不能求解' return endendx = zeros(n, 1);

18、 k = 0; kend = 10000; r = 1; while k< = kend&r>e x0 = x; for i = 1:n s = 0; for j = 1:i - 1s = s + H(i, j) * x0(j); end for j = i + 1:ns = s + H(i, j) * x0(j); end x(i) = b(i) / H(i, i) - s / H(i, i); end r = norm(x - x0, inf); k = k + 1; end if k>kend'迭代不收敛，失败' else'求解成功

19、9;x k end输出结果为：ans = 迭代不收敛，失败6.2 用SOR迭代法求解建立新的m文件unit62.m：代码如下：function s = unit62 (n, w); H = hilb(n); b = H * ones(n, 1); e = 0.00001; for i = 1:n if H(i, i) = = 0 '对角元为零，不能求解' return endendx = zeros(n, 1); k = 0; kend = 10000; r = 1; while k< = kend&r>e x0 = x; for i = 1:n s = 0

20、; for j = 1:i - 1s = s + H(i, j) * x(j); end for j = i + 1:ns = s + H(i, j) * x0(j); end x(i) = (1 - w) * x0(i) + w / H(i, i) * (b(i) - s); end r = norm(x - x0, inf); k = k + 1; end if k>kend'迭代不收敛，失败' else '求解成功'xk end输入：unit62(6, 1)输出结果为：ans = 求解成功x = 0.99931.01310.95371.03741.0

21、2960.9662k = 997ans = 0.6487当n = 6，w = 1.25时，ans = 求解成功，x = 0.9992, 1.0131, 0.9558, 1.0375, 1.0197, 0.9741，k = 3121，ans = 0.6480；当n = 6，w = 1.5时，ans = 求解成功，x = 0.9995, 1.0081, 0.9726, 1.0232, 1.0123, 0.9839，k = 5635，ans = 0.6471；当n = 8，w = 1时，ans = 求解成功，x = 0.9996, 1.0039, 0.9938, 0.9824, 1.0071，1.0

22、200，1.0093, 0.9790 ，k = 3808，ans = 0.6601；当n = 8，w = 1.25时，ans = 求解成功x = 0.9998, 1.0004, 1.0103, 0.9685, 1.0109, 1.0201, 1.0100, 0.9795，k = 3949，ans = 0.6601；当n = 8，w = 1.5时，ans = 求解成功，x = 1.0000, 0.9975, 1.0231, 0.9469, 1.0275, 1.0118, 1.0158, 0.9771，k = 4031，ans = 0.6602；当n = 10，w = 1时，ans = 求解成功，

23、x = 1.0001, 0.9952, 1.0275, 0.9677, 0.9816, 1.0093, 1.0241, 1.0208, 1.0018, 0.9713，k = 2673，ans = 0.6677, ；当n = 10，w = 1.25时，ans = 求解成功，x = 1.0004, 0.9904, 1.0466, 0.9421, 0.9886, 1.0123, 1.0272, 1.0214, 1.0010, 0.9694，k = 2562，ans = 0.6677；当n = 10，w = 1.5时，ans = 求解成功，x = 1.0006, 0.9857, 1.0688, 0.9

24、028, 1.0177, 0.9993, 1.0377, 1.0169, 1.0030, 0.9669 ，k = 2463，ans = 0.66797 用MATLAB求解非线性方程的根题目求下列方程的实根(1) ; (2) .要求：(1) 设计一种不动点迭代法，要使迭代序列收敛，然后再用斯特芬森加速迭代，计算到为止.(2) 用牛顿迭代，同样计算到.输出迭代初值及各次迭代值和迭代次数k，比较方法的优劣.7.1 不动点迭代法建立新的m文件unit71.m：代码如下：function y, n = unit71(x, eps)if nargin = = 1 eps = 1.0e - 8; else

25、if nargin<1 error returnendx1 = gg(x); n = 1; while (norm(x1 - x)> = 1e - 8)&&(n< = 10000) x = x1; x1 = gg(x); n = n + 1; endy = x; 若是function f = gg(x)f(1) = (exp(x) - 2 - x2) / ( - 3); 在在命令窗口输入：unit71(1)输出结果为：n = 15ans = 0.2575若是function f = gg(x)f(1) = - (x3 + 2 * x2 - 20) / 10;

26、在命令窗口输入：unit71(1)输出结果为：n = 10001ans = 0.54897.2 斯特芬森加速迭代法建立新的m文件unit72.m：代码如下：function y = unit72(x0, eps)if nargin<2 eps = 1.0e - 8; enda = g(x0); b = g(a); x1 = x0 - (a - x0)2 / (b - 2 * a + x0); n = 1; while (abs(x1 - x0)> = eps)&(n< = 10000) x0 = x1; a = g(x0); b = g(a); x1 = x0 - (a - x0)2 / (b - 2 * a + x0); n = n + 1; endx1n若是function y = g(x); y = (exp(x) - 2 - x2) / ( - 3); 在命令窗口输入：unit72(1)输出结果为：x1 = 0.2575n = 4若是function y = g(x); y = - (x3 + 2 * x2 - 20) / 10在命令窗口输入：unit72(1)输出结果为：x1 = 1.3688n = 57.3 牛顿迭代法建立新的m文件unit73.m：代码如下：function y = unit73(x0, eps)if narg