数值分析法求正弦余弦积分函数

天津职业技术师范大学天津职业技术师范大学 课课 程程 设设 计计 任任 务务 书书 理 学院 数学 1403 班 学生 张群 课程设计课题 用数值积分法计算正弦积分函数和余弦积分函数 一 课程设计工作日自 2016 年 7 月 4 日至 2016 年 7 月 5 日 二 同组学生 无 三 课程设计任务要求 包括课题来源 类型 目的和意义 基本要求 完成 时间 主要参考资料等 课题来源 课题来源 教师自拟 类型 类型 理论研究 目的和意义 目的和意义 培养学生对数值分析中主要算法的应用能力 探索相关算法 之间的内在联系 基本要求 基本要求 根据数值分析课程所学的知识 应用 Matlab 软件编写程序 完 成对算法及其内在原理的实验研究 完成时间 完成时间 参考资料 参考资料 数值分析 李庆扬等 清华大学出版社 指导教师签字 指导教师签字 教研室主任签字 教研室主任签字 一 问题叙述 用数值积分法计算正弦积分函数和余弦积分函数用数值积分法计算正弦积分函数和余弦积分函数 提示 正弦积分 余弦函数 0 sin x t si xdt t cos x t ci xdt t 要求 1 编写函数 对任意给定的 x 可求出值 2 使用尽可能少的正余弦函数的调用 计算更精确的值 用多种方法 创新方法 2 问题分析 1 数值积分基本原理 用数值分析求解积分的数值方法有很多 如简单的梯形法 矩形法 辛普森 Simpson 法 牛顿 科斯特 Newton Cotes 法等都是常用的方法 它们的基本思想都是将整 个积分区间 a b 分成 n 个子区间 xi xi 1 i 1 2 n 其中 x1 a xn 1 b 这样求定积分问题就分解为求和问题 2 本题要求用数值积分法计算正弦积分函数和余弦函数积分 那么 应该从编写函数的入手 建立 function 的 m 文件 通过对函数的调 用 对任意跟定的 x 值 求出积分函数的值 3 数值积分法求解问题 1 1 梯形公式 矩形公式梯形公式 矩形公式 首先 积分中值定理告诉我们 在积分区间 a b 内存在一点 成 立dx b a f 就是说 底为 b a 而高为 f 的矩 b a xf 形面积恰等于所求区边梯形的面积 如果我们用两端点 高度 f a 与 f b 的算术平均值作为平均高度 f 的近似值 这样 导出的求积公式dx f a f b 便是我们熟悉的梯形 b a xf 2 a b 公式 将积分区间 a b n 等分 每个小区间宽度均为 h h 称 n a b 为积布步长 记 a x0 x1 xk xo b 在小区间 上 用 小 矩 形 面 积 近 似 小 曲 边 梯 形 的 面 积 若 分 别 取 左 端 点 和 右 端 点 的 函 数 值 为 小 矩 形 的 高 则 分 别得到两个曲边 梯形面积的近似计算公式 具体程序如下 clear x linspace 0 pi dx x 2 y sin x s1 sum y dx s2 trapz y dx sc1 cumsum y dx sc2 cumtrapz y dx plot x cos x 1 x sc1 x sc2 o hold on 00 511 522 533 5 0 0 2 0 4 0 6 0 8 1 1 2 1 4 1 6 1 8 2 由图可知这种方法精度太低 应选择其他方法 2 2 quadquad 函数 函数 quan1quan1 函数函数 正弦 function y si t a 1e 8 函数在 0 点无界 去掉 0 点 y quad sin x x a t y quadl sin x x a t 余弦 function y ci t a 1e1 函数在 0 点无界 去掉 0 点 y quad cos x x a t y quadl cos x x a t 图像 x 1 100 for i 1 100 y2 i si x i end plot x y2 r title 辛普森 0102030405060708090100 0 9 1 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 化 化 化 x 1 100 for i 1 100 y2 i ci x i end plot x y2 b title 辛普森 0102030405060708090100 400 200 0 200 400 600 800 1000 1200 1400 化 化 化 给定任意 x 值 均可计算出对应的正弦 余弦函数积分 但从结 果可以看出精度不是很高 3 复合求积公式复合求积公式 由于牛顿 科特斯公式在 n 8 时不具有稳定性 故不可能通过提 高阶的方法来提高求积精度 为了提高精度通常可把积分区间分成 若干子区间 通常是等分 再在每个子区间上用低级求积公式 这种方法为复合求积法 3 3 1 复合梯形公式 将区间划分为 n 等分 分点在每个 ba 1 0 nk n ab hkhaxk 子区间上采用梯形公式 则得 1 1 0 1 nkxx kk 2 1 1 0 1 0 1 fRxfxf h dxxfdxxfI nk n k k b a n k x x k k 记 1 1 1 1 0 2 22 n k bkk n k kn xfxfaf h xfxf h T 称为复合梯形公式 复合梯形公式 复合梯形公式的余项 1 1 0 3 12 kkk n k knn xxf h TIfR 由于且 2 baCxf max 1 min 10 1 0 10 nk k n k k nk f n f 所以使 ba k n k f n f 1 0 1 于是复合梯形公式的余项为 2 12 fh ab fRn 事实上只要设 则可得收敛性 只要把改写成为 baCxf n T n k k n k kn xf n ab xf n ab T 1 1 0 2 1 程序如下 正弦 function T n fhtxs a b n h b a n for k 0 n x k 1 a k h if x k 1 0 x k 1 10 10 end end T 1 h 2 SS x 1 SS x n 1 for i 2 n F i h SS x i end T 2 sum F T n T 1 T 2 余弦 function T n fhtxc a b n h b a n for k 0 n x k 1 a k h if x k 1 0 x k 1 10 10 end end T 1 h 2 CC x 1 CC x n 1 for i 2 n F i h CC x i end T 2 sum F T n T 1 T 2 图像 正弦 余弦 01020304050 0 8 1 1 2 1 4 1 6 1 8 2 化 化 化 化 020406080100 2 0 2 4 6 8 10 x 10 10 3 3 2 复合新普斯求积公式 将区间划分为 n 等分 在每个子区间上采用辛普森公式 若记 ba 1 kk xx 则得 2 1 21 hxx kk 1 0 n k b a dxxfdxxfI 4 6 1 0 21 fRxfxfxf h n n k kkk 称为复合辛普森求积公式 复合辛普森求积公式 程序如下 正弦 function S n fhxpss a b n h b a n for k 0 n x k 1 a k h x k k 1 x k 1 1 2 h if x k 1 0 x k k 1 0 x k 1 10 10 x k k 1 10 10 end end S 1 h 6 SS x 1 SS x n 1 for i 2 n F 1 i h 3 SS x i end for j 1 n F 2 j 2 h 3 SS x k j end S 2 sum F 1 sum F 2 S n S 1 S 2 余弦 function S n fhxpsc a b n h b a n for k 0 n x k 1 a k h x k k 1 x k 1 1 2 h if x k 1 0 x k k 1 0 x k 1 10 10 x k k 1 10 10 end end S 1 h 6 CC x 1 CC x n 1 for i 2 n F 1 i h 3 CC x i end for j 1 n F 2 j 2 h 3 CC x k j end S 2 sum F 1 sum F 2 S n S 1 S 2 图像与复合梯形所得图像基本相同 深入分析两只复合函数的优劣 对于积分函数 假设 x 1 则将区间 0 1 划分为 8 等份 0 sin x t si xdt t 应用复合梯形求得 T8 0 9456909 而如果将 0 1 分为 4 等份 应用复合辛普森有 S4 0 9460832 通过参考数值分析 李庆阳 的结论 发现无论是复合梯形公式还是复合辛普 森公式 最终结果都会随着 h 值的减小而更加精确 对复合梯形公式和复合辛 普森公式计算出的结果进行比较 发现复合梯形法的结果T8只有两位有效数 字 而复合辛普森的结果却有六位有效数字 所以复合辛普森公式计算出的结 果更加的精确 4 插值型的求积公式插值型的求积公式 clc clear x0 0 0 5 5 y0 Inf 1 7552 0 5403 0 0472 0 2081 0 3205 0 3300 0 2676 0 1634 0 0468 0 0567 所求积分函数的数值 pp csape x0 y0 默认的边界条件 Lagrange 边界条件 format long g chazhi pp coefs 显示每个区间上三次多项式的系数 s quadl t ppval pp t 0 5 求积分 format 恢复短小数的显示格式 x 0 0 1 5 y cos x x y1 spline x0 y0 x z 0 x hold on plot x z x y k x y1 r plot x0 y0 hold off clear x0 0 0 5 5 y0 Inf 1 7552 0 5403 0 0472 0 2081 0 3205 0 3300 0 2676 0 1634 0 0468 0 0567 所求积分函数的数值 pp csape x0 y0 默认的边界条件 Lagrange 边界条件 format long g chazhi pp coefs 显示每个区间上三次多项式的系数 s quadl t ppval pp t 0 5 求积分 format 恢复短小数的显示格式 x 0 0 1 5 y cos x x y1 spline x0 y0 x z 0 x hold on plot x z x y k x y1 r plot x0 y0 hold off 如图所示 5 高斯求积公式 function ql Ak xk gsqj fun a b n tol if nargin 1 a 1 b 1 n 7 tol 1e 8 elseif nargin 3 n 7 tol 1e 8 elseif nargin 4 tol 1e 8 elseif nargin 2 nargin 5 error The Number of Input Arguments Is Wrong end 计算求积节点 syms x p sym2poly diff x 2 1 n 1 n 1 2 n factorial n tk roots p 求积节点 计算求积系数 Ak zeros n 1 1 for i 1 n 1 xkt tk xkt i pn poly xkt fp x polyval pn x polyval pn tk i Ak i quadl fp 1 1 tol 求积系数 end 积分变量代换 将 a b 变换到 1 1 xk b a 2 tk b a 2 检验积分函数 fun 有效性 fun fcnchk fun vectorize 计算变量代换之后积分函数的值 SS fu