《数值计算课程7090307龙凯翔.doc》由会员分享，可在线阅读，更多相关《数值计算课程7090307龙凯翔.doc（16页珍藏版）》请在金锄头文库上搜索。
1、东北大学秦皇岛分校数值计算课程设计报告系 别信息与计算科学专 业学 号7090307姓 名龙凯翔指导教师张建波 姜玉山成 绩教师评语：指导教师签字： 2011年7月13日 信息与计算科学系数值计算课程设计报告 第 3 页1 绪 论数值分析是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。MATLAB是主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计
2、算的众多科学领域提供了一种全面的解决方案。2 用MATLAB进行插值计算题目下列数据点的插值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 =
3、polyval(p, xi); subplot(2, 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(插值点, 三次样条插值)xi = 0:0.01:1; p = polyfit(x, y, 8); yi = polyval(p, xi); subplot(2, 2, 3); plot(x
4、i, 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(三次样条插值) 图12.2 图样运行文件得到图12.3 分析比较两个函数的图像可知，在区间0, 64上三次样条插值函数要更精确一些，在区间0, 1上拉格朗日插值函数仍然不如三次样条插值函数更精确。3 多项式曲线拟合题目由实验给出数据表X0.00.10.20.30.50.81.0Y1.00.410.500.610.
5、912.022.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
6、(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轴); legend(插值点, 4次曲线拟合)poly2str(a, x)poly2str(b, x) 3.2 图形和结果运行程序结果如下：ans = - 6.6221 x3 + 12.8147 x2 - 4.6591 x + 0.92659ans = 2.8853 x4 - 12.3348 x3 + 16.2747
7、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(x轴); ylabel(y轴); legend(插值点, 5次曲线拟合)poly2str(a, x)运行程序结果如
8、下：ans = - 79.3261 x5 + 195.4554 x4 - 172.7104 x3 + 69.0498 x2 - 11.0044 x + 0.99547图像如下：图3由图像可知，次数越高，拟合越准确。4 数值积分题目用不同数值方法计算积分.(1)取不同的步长h.分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能呢个再被改善？(2)用龙贝格求积计算完成问题(1).(3)用自适应辛普森积分，使其精度达到4.1 复合梯形公式代码建立新的m文件unit41.m：代码如下：clc, clear; x = 0.
9、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.4875e - 005；h = 0.0001，Fx = - 0.4444，error = - 2.0338e - 006；h = 0.00001，Fx = - 0.4444，error = - 6.4496e - 0
10、08.可见，没有一个最小的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 = 0.00001，Fx = - 0.4444，error = - 6.3537e - 005；h = 0.000001，Fx = - 0.4444，error = - 2.9938e - 006；h = 0.0000001，Fx
11、= - 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之间。4.3龙贝格算法建立新的m文件unit43.m：代码如下：function q, n = unit43(f, a, b); M = 1; abs0 = 10; k = 0; T =
12、 zeros(1, 1); h = b - a; T(1, 1) = (h / 2) * (subs(f, a) + subs(f, b); while abs00.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) = (4j) * T(k + 1, j) - T(k, j) / (4j - 1); end a
13、bs0 = 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(inline(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 10det(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.30
