数值计算课程7090307龙凯翔

数值计算课程7090307龙凯翔 东北大学秦皇岛分校 数值计算课程 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 报告 软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载 系 别 信息与计算科学 专 业 数学与应用数学 学 号 7090307 姓 名 龙凯翔 指导教师 张建波 姜玉山 成 绩 教师评语: 指导教师签字: 2011年7月13日 信息与计算科学系数值计算课程设计报告 第 1 页 1 绪 论 数值分析是研究分析用计算机求解数学计算问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。 MATLAB是主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决 方案 气瓶 现场处置方案 .pdf气瓶 现场处置方案 .doc见习基地管理方案.doc关于群访事件的化解方案建筑工地扬尘治理专项方案下载 。 2 用MATLAB进行插值计算 题目 下列数据点的插值 X 0 1 4 9 16 25 36 49 64 Y 0 1 2 3 4 5 6 7 8 可以得到平方根函数的近似，在区间[0, 64]上作图。 Lx()8(1) 用这9个点作8次多项式插值. Sx()(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, 2, 1); plot(x, y, 'x', xi, yi, 'k') xlabel('x轴'); ylabel('y轴'); legend('插值点', '8次差值多项式') 信息与计算科学系数值计算课程设计报告 第 2 页 yi = interp1(x, y, xi, 'cubic'); subplot(2, 2, 2); plot(x, y, 'x', xi, yi, 'k') xlabel('x轴'); ylabel('y轴'); legend('插值点', '三次样条插值') 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轴'); ylabel('y轴'); legend('三次样条插值') 图1 信息与计算科学系数值计算课程设计报告 第 3 页 2.2 图样 运行文件得到图1 2.3 分析 比较两个函数的图像可知，在区间[0, 64]上三次样条插值函数要更精确一些，在区间[0, 1]上拉格朗日插值函数仍然不如三次样条插值函数更精确。 3 多项式曲线拟合 题目 由实验给出数据 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 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 试求3次、4次多项式的曲线拟合，再根据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。 3.1程序文件 建立新的m文件unit3.m: 代码如下: 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, 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('插值点', '3次曲线拟合') subplot(1, 2, 2); plot(x, y, 'x', xi, bb, 'k') xlabel('x轴'); ylabel('y轴'); 信息与计算科学系数值计算课程设计报告 第 4 页 legend('插值点', '4次曲线拟合') poly2str(a, 'x') poly2str(b, 'x') 3.2 图形和结果 运行程序结果如下: ans = - 6.6221 x^3 + 12.8147 x^2 - 4.6591 x + 0.92659 ans = 2.8853 x^4 - 12.3348 x^3 + 16.2747 x^2 - 5.2987 x + 0.94272 图像如下: 图2 3.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('x轴'); ylabel('y轴'); legend('插值点', '5次曲线拟合') poly2str(a, 'x')运行程序结果如下: 信息与计算科学系数值计算课程设计报告 第 5 页 ans = - 79.3261 x^5 + 195.4554 x^4 - 172.7104 x^3 + 69.0498 x^2 - 11.0044 x + 0.99547 图像如下: 图3 由图像可知，次数越高，拟合越准确。 4 数值积分 题目 14lnxxdx,,,09用不同数值方法计算积分. (1)取不同的步长h.分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能呢个再被改善, (2)用龙贝格求积计算完成问题(1). ,410(3)用自适应辛普森积分，使其精度达到 4.1 复合梯形公式代码 建立新的m文件unit41.m: 信息与计算科学系数值计算课程设计报告 第 6 页 代码如下: clc, clear; x = 0.0000001:0.1:1; y = x.^(1 / 2). * log(x); Fx = trapz(x, y) error = Fx + 4 / 9 运行文件结果为: Fx = - 0.4123 error = 0.0321 取不同步长h: h = 0.01，Fx = - 0.4431，error = - 0.0014; h = 0.001，Fx = - 0.4444，error = - 5.4875e - 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.4440 error = 4.3273e - 004 取不同的步长h: h = 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.0000000000001与 0.00000000000001之间。 信息与计算科学系数值计算课程设计报告 第 7 页 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(k + 1, 1) = T(k, 1) / 2 + h * p; M = 2 * M; for j = 1:k T(k + 1, j + 1) = ((4^j) * T(k + 1, j) - T(k, j)) / (4^j - 1); end abs0 = abs(T(k + 1, j + 1) - T(k, j)); end q = T(k + 1, k + 1); n = k; 在命令窗口输入: [Fx, n] = unit43('sqrt(x) * log(x)', 10^( - 8), 1) 输出结果为: Fx = - 0.4444 n = 9 4.4 自适应辛普森算法 信息与计算科学系数值计算课程设计报告 第 8 页 Fx = quad(inline('x.^(1 / 2). * log(x)'), 0.000001, 1) 输出结果为: Fx = - 0.4444。 5 解线性方程 题目 线性方程组Ax = b的A及b为 1078732,,,,,,,,756523,,,,A = ，b = ， ,,,,8610933,,,,,,,,7591031,,,, ，，1,1,1,12则解x = .用MATLAB内部函数求detA及A的所有特征值和cond(A) . 若令 1078.17.2,,,,7.085.0465,,A + δA = ， ,,85.989.899,,,,6.99599.98,, δx2求解(A + δA)(x + δx) = b, 输出向量δx和.从理论结果和实际计算两方面分 δx/xδA/A2222析线性方程组Ax = b解的相对误差及A的相对误差的关系. 25.1 求detA及A的所有特征值和cond(A) 输入程序: A = [10 7 8 7; 7 5 6 5; 8 6 10 9; 7 5 9 10] det(A) [V, D] = eig(A) norm(A, 2) 输出结果为: A = 10 7 8 7 7 5 6 5 8 6 10 9 7 5 9 10 信息与计算科学系数值计算课程设计报告 第 9 页 ans = 1 V = 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.5209 D = 0.0102 0 0 0 0 0.8431 0 0 0 0 3.8581 0 0 0 0 30.2887 ans = 30.2887 则得到detA = 1 A的特征值为0.0101 0.8431 3.8581 30.2887 condA，，2 = 30.2887 6 解线性方程组的迭代法 题目 Hx,bHnn给出线性方程组其中系数矩阵为希尔伯矩阵: 1h,n，nijH,(h),Ri，j,1nij，，i = 1, 2……n. **Tn，nb,Hxx,(1,1？,1),Rn假设，，若取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 信息与计算科学系数值计算课程设计报告 第 10 页 if H(i, i) = = 0 '对角元为零，不能求解' return end end x = zeros(n, 1); k = 0; kend = 10000; r = 1; while k< = kend&r>e x0 = x; for i = 1:n s = 0; for j = 1:i - 1 s = s + H(i, j) * x0(j); end for j = i + 1:n s = 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 '求解成功' x k end 输出结果为: ans = 迭代不收敛，失败 6.2 用SOR迭代法求解 建立新的m文件unit62.m: 代码如下: function s = unit62 (n, w); H = hilb(n); 信息与计算科学系数值计算课程设计报告 第 11 页 b = H * ones(n, 1); e = 0.00001; for i = 1:n if H(i, i) = = 0 '对角元为零，不能求解' return end end x = zeros(n, 1); k = 0; kend = 10000; r = 1; while k< = kend&r>e x0 = x; for i = 1:n s = 0; for j = 1:i - 1 s = s + H(i, j) * x(j); end for j = i + 1:n s = 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 '求解成功' x k end 输入:unit62(6, 1) 输出结果为: ans = 求解成功 x = 0.9993 1.0131 0.9537 1.0374 1.0296 0.9662 信息与计算科学系数值计算课程设计报告 第 12 页 k = 997 ans = 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.0200， 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 = 求解成功，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.9028, 1.0177, 0.9993, 1.0377, 1.0169, 1.0030, 0.9669] ，k = 2463，ans = 0.6679 7 用MATLAB求解非线性方程的根 题目 求下列方程的实根 2xx,3x，2,e,0(1) ; 32x，2x，10x,20,0(2) . 要求:(1) 设计一种不动点迭代法，要使迭代序列收敛，然后再用斯特芬森加速迭代， ,8,8x,x,10x,x,10kkkk,1,1计算到为止.(2) 用牛顿迭代，同样计算到.输出迭代初值及 各次迭代值和迭代次数k，比较方法的优劣. 7.1 不动点迭代法 建立新的m文件unit71.m: 信息与计算科学系数值计算课程设计报告 第 13 页 代码如下: function [y, n] = unit71(x, eps) if nargin = = 1 eps = 1.0e - 8; else if nargin<1 error return end x1 = gg(x); n = 1; while (norm(x1 - x)> = 1e - 8)&&(n< = 10000) x = x1; x1 = gg(x); n = n + 1; end y = x; 若是function f = gg(x) f(1) = (exp(x) - 2 - x^2) / ( - 3); 在在命令窗口输入:unit71(1) 输出结果为: n = 15 ans = 0.2575 若是function f = gg(x) f(1) = - (x^3 + 2 * x^2 - 20) / 10; 在命令窗口输入:unit71(1) 输出结果为: n = 10001 ans = 0.5489 7.2 斯特芬森加速迭代法 建立新的m文件unit72.m: 代码如下: function y = unit72(x0, eps) if nargin<2 eps = 1.0e - 8; end a = g(x0); b = g(a); x1 = x0 - (a - x0)^2 / (b - 2 * a + x0); n = 1; 信息与计算科学系数值计算课程设计报告 第 14 页 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; end x1 n 若是function y = g(x); y = (exp(x) - 2 - x^2) / ( - 3); 在命令窗口输入:unit72(1) 输出结果为: x1 = 0.2575 n = 4 若是function y = g(x); y = - (x^3 + 2 * x^2 - 20) / 10 在命令窗口输入:unit72(1) 输出结果为: x1 = 1.3688 n = 5 7.3 牛顿迭代法 建立新的m文件unit73.m: 代码如下: function y = unit73(x0, eps) if nargin<2 eps = 1.0e - 8; end x1 = x0 - fa(x0) / faa(x0); n = 1; while (abs(x1 - x0)> = eps)&(n< = 10000) x0 = x1; x1 = x0 - fa(x0) / faa(x0); n = n + 1; end x1 n 若是function y = fa(x); y = x^2 - 3 + 2 * x - exp(x); function y = faa(x); y = 2 * x - 3 - exp(x); 在命令窗口输入:unit73(0) 信息与计算科学系数值计算课程设计报告 第 15 页 输出结果为: x1 = - 3.0123 n = 33 若是function y = fa(x); y = x^3 + 2 * x^2 + 10 * x - 20; function y = faa(x); y = 3 * x^2 + 4 * x + 10; 在命令窗口输入:unit73(0) 输出结果为: x1 = 1.3688 n = 6 通过三种方法的比较，可以看出斯特芬森加速迭代法迭代次数最少，方法最优，牛顿迭代法次之，而不动点迭代法次数最多，不动点迭代法针对收敛函数求解，所以在不动点迭代法求解之前要先判断函数敛散性。牛顿迭代法因其效率高而应用最为广泛，但是它对重很收敛很慢，而且对初值要求严格。 结 论 数值分析是研究如何用计算机解决实际问题的课程(将MATLAB与数值分析课程结合起来，开阔了思路，拓展了解决问题的方法，取得了较好的效果(MATLAB在数值计算中发挥着十分重要的作用( 通过这段时间对MATLAB的学习，我对MATLAB有了更进一步的了解，能够运用MATLAB来编写简单的函数和程序，也能进行符号计算，尤其是我可以通过MATLAB来解决数值分析的一些问题(我对计算机解决实际问题有了一个全新的认识(我会继续努力地学习MATLAB软件，掌握更多的相关知识( 参考文献 [1] 李庆杨, 王能超, 易大义. 数值分析(第5版)[M]. 北京: 清华大学出版社, 2006. [2] 胡良剑, 孙晓君, MATLAB数学实验 [M]. 北京: 高等教育出版社,2006.