matlab实现数值分析插值及积分

Matlab 实现数值分析插值及积分实现数值分析插值及积分 摘要 摘要 数值分析 numerical analysis 是研究分析用计算机求解数学计算问题的数值计算方法及 其理论的学科 是数学的一个分支 它以数字计算机求解数学问题的理论和方法为研究对 象 在实际生产实践中 常常将实际问题转化为数学模型来解决 这个过程就是数学建模 学习数值分析这门课程可以让我们学到很多的数学建模方法 分别运用matlab数学软件编程来解决插值问题和数值积分问题 题目中的要求是计算 差值和积分 对于问题一 可以分别利用朗格朗日插值公式 牛顿插值公式 埃特金逐次 线性插值公式来进行编程求解 具体matlab代码见正文 编程求解出来的结果为 其中Aitken插值计算的结果图如下 对于问题二 可以分别利用复化梯形公式 复化的辛卜生公式 复化的柯特斯公式编 写程序来进行求解 具体matlab代码见正文 编程求解出来的结果为 0 6932 其中复化梯形公式计算的结果图如下 问题重述问题重述 问题一 问题一 已知列表函数 表格 1 01234 121782257 分别用拉格朗日 牛顿 埃特金插值方法计算 问题二 问题二 用复化的梯形公式 复化的辛卜生公式 复化的柯特斯公式计算积分 使精度小 于 5 问题解决问题解决 问题一 插值方法问题一 插值方法 对于问题一 用三种差值方法 拉格朗日 牛顿 埃特金差值方法来解决 一 拉格朗日插值法 一 拉格朗日插值法 拉格朗日插值多项式如下 首先构造个插值节点上的插值基函数 对任一点所对应的插值基函1 n n xxx 10 n i x 数 由于在所有取零值 因此有因子 xli 1 1 1 0 niijxj xli 又因是一个次数不超过的多项式 所以 110nii xxxxxxxx xlin 只可能相差一个常数因子 固可表示成 xli 110niii xxxxxxxxAxl 利用得 1 ii xl 1 110niiiiii xxxxxxxx A 于是 2 1 0 110 110 ni xxxxxxxx xxxxxxxx xl niiiiii nii i 因此满足 的插值多项式可表示为 iin yxL 2 1 0 ni n j jjn xlyxL 0 从而次拉格朗日插值多项式为 n 2 1 0 0 nixlyxL n j ijjin matlab 编程 编程 编程思想 编程思想 主要从上述朗格朗日公式入手 依靠循环 运用 poly 函数和 conv 函数 表示拉格朗日公式 其中的 poly i 函数表示以 i 作为根的多项式的系数 例如 poly 1 表示 x 1 的系数 输出为 1 1 而 poly poly 1 表示 x 1 x 1 x 2 2 x 1 的系 数 输出为 1 2 1 而 conv 表示多项式系数乘积的结果 例如 conv poly 1 poly 1 输出为 1 2 1 所以程序最后结果为 x n x n 1 x 2 x 1 n 的值据结果 的长度为准 的对应系数 在命令窗口输入 edit lagran 来建立 lagran m 文件 文件中的程序如下 function c l lagran x y w length x n w 1 l zeros w w for k 1 n 1 v 1 for j 1 n 1 if k j v conv v poly x j x k x j end end l k v end c y l 输入 输入 x 0 1 2 3 4 y 1 2 17 82 257 lagran x y 运行结果为运行结果为 ans 1 0000 0 0000 0 0000 0 1 0000 结果为 结果为 如图表如图表1 1 图表 1 二 牛顿插值法二 牛顿插值法 newton 插值多项式的表达式如下 010011 nnn Nxccxxcxxxxxx 其中每一项的系数 ci 的表达式如下 12011 01 0 ii ii i f x xxf x xx cf x xx xx 即为 f x 在点处的 i 阶差商 由差商 01 i x xx ii f xf x 1 2 in 的性质可知 01 i f x xx 01 00 1 ii ij jk jk kj f x xxf x xx matlab 编程 编程 编程思想 编程思想 主要从上述牛顿插值公式入手 依靠循环 运用 poly 函数和 conv 函数 表示拉格朗日公式 其中的 poly i 函数表示以 i 作为根的多项式的系数 例如 poly 1 表示 x 1 的系数 输出为 1 1 而 poly poly 1 表示 x 1 x 1 x 2 2 x 1 的系 数 输出为 1 2 1 而 conv 表示多项式系数乘积的结果 例如 conv poly 1 poly 1 输出为 1 2 1 所以程序最后结果为 x n x n 1 x 2 x 1 n 的值据结果 的长度为准 的对应系数 在命令窗口输入 edit nowpoly 来建立 newpoly m 文件 文件中的程序如下 function c d newpoly x y n length x d zeros n n d 1 y for j 2 n for k j n d k j d k j 1 d k 1 j 1 x k x k j 1 end end c d n n for k n 1 1 1 c conv c poly x k m length c c m c m d k k end 输入 输入 x 0 1 2 3 4 y 1 2 17 82 257 newpoly x y 运行结果为运行结果为 ans 1 0 0 0 1 所以 如图表如图表2 2 图表 2 三 埃特金插值法 三 埃特金插值法 Aitken 插值公式如下 递推表达式为 当 n 1 时 当 n 2 时 其中的带入递推表达式求得 由此递推下去 最终得到的结果 matlab 编程 编程 编程思想 编程思想 埃特金插值多项式又称作 Aitken 逐次线性插值多项式 根据公式的特点 可 以利用 2 次嵌套循环将公式表示出来 在命令窗口输入 edit Aitken 来建立 Aitken m 文件 文件中的程序如下 function f Aitken x y syms z n length x y1 1 n z for i 1 n 1 for j i 1 n y1 j y j z x i x j x i y i z x j x i x j end y y1 simplify y1 end simplify y1 n f collect y1 n 输入 输入 x 0 1 2 3 4 y 1 2 17 82 257 Aitken x y 运行结果为运行结果为 ans z 4 1 所以 如图表如图表 3 图表 3 问题二 复化积分问题二 复化积分 对于问题二来说 用复化的梯形公式 复化的辛卜生公式 复化的柯特斯公式结局问 题 计算积分 使精度小于 5 一一 复化的梯形公式 复化的梯形公式 复化梯形的迭代公式为 matlab 编程 编程 程序程序 1 求 f x 的 n 阶导数 在命令窗口输入 edit qiudao 来建立 qiudao m 文件 文件中的程序如下 function d qiudao x n syms x f 1 x n input 输入导数阶数 d diff f x n 输入输入 qiudao x n 输入所求导数阶数 2 显示 显示 n 2 ans 2 x 3 结果为结果为 f2 2 x 3 如图表4 图表 4 程序程序 2 在命令窗口输入 edit tixing 来建立 tixing m 文件 文件中的程序如下 function y tixing syms x 定义自变量定义自变量 x f inline 1 x x 定义函数定义函数 f x 1 x f2 inline 2 x 3 x 定义定义 f x 的二阶导数的二阶导数 输入程序输入程序 1 里求出的里求出的 f2 即可 即可 f3 2 x 3 因因 fminbnd 函数求的是表达式的最小值 且要求表 函数求的是表达式的最小值 且要求表 达式带引号 故取负号 以便求最大值达式带引号 故取负号 以便求最大值 e 5 10 5 精度要求值精度要求值 a 1 积分下限积分下限 b 2 积分上限积分上限 x1 fminbnd f3 1 2 求负的二阶导数的最小值点 也就是求二阶导数的最大求负的二阶导数的最小值点 也就是求二阶导数的最大 值点对应的值点对应的 x 值值 for n 2 1000000 求等分数求等分数 n Rn b a 12 b a n 2 f2 x1 计算余项计算余项 if abs Rn e 用余项进行判断用余项进行判断 break 符合要求时结束符合要求时结束 end end h b a n 求求 h Tn1 0 for k 1 n 1 求连加和求连加和 xk a k h Tn1 Tn1 f xk end Tn h 2 f a 2 Tn1 f b fprintf 用复化梯形算法计算的结果 Tn disp Tn fprintf 等分数 n disp n 输出等分数输出等分数 输入 输入 tixing 运行结果为运行结果为 用复化梯形计算的结果 Tn 0 6932 等分数 n 58 结果如图表结果如图表 5 图表 5 二二 复化的复化的辛卜生辛卜生公式 公式 复化 simpson 迭代公式为 matlab 编程 编程 程序程序 1 求 f x 的 n 阶导数 在命令窗口输入 edit qiudao 来建立 qiudao m 文件 文件中的程序如下 function d qiudao x n syms x f 1 x n input 输入导数阶数 d diff f x n 输入输入 qiudao x n 输入所求导数阶数 4 显示 显示 n 4 ans 24 x 5 结果为结果为 f2 24 x 5 如图表6 图表 6 程序程序 2 在命令窗口输入 edit xinpusheng 来建立 xinpusheng m 文件 文件中的程序如下 function y xinpusheng syms x f inline 1 x x f2 inline 24 x 5 x f3 24 x 5 e 5 10 5 a 1 b 2 x1 fminbnd f3 1 2 for n 2 1000000 Rn b a 180 b a 2 n 4 f2 x1 if abs Rn clear xinpusheng 运行结果为运行结果为 用Simpson公式计算的结果 Sn 0 6932 等分数 n 4 结果如图表结果如图表 7 图表 7 三三 复化的复化的柯特斯柯特斯公式 公式 牛顿 柯特斯公式如下 matlab 编程 编程 1 在命令窗口输入 edit NewtonCotes 来建立 NewtonCotes m 文件 文件中的程序如下 function y Ck Ak NewtonCotes fun a b n if nargin 1 mm nn size fun if mm 8 error 为了保证 NewtonCotes 积分的稳定性 最多只能有 9 个等距节点 elseif nn 2 error fun 构成应为 第一列为 x 第二列为 y 并且个数为小于 10 的等距节点 end 0 n b kk a k fx dxA fx 011 011 bb kkn kk aa kkkkkkn xxxxxxxx Alx dxdx xxxxxxxx xk fun 1 fk fun 2 a min xk b max xk n mm 1 elseif nargin 4 xk linspace a b n 1 if isa fun function handle fx fun xk else error fun 积分函数的句柄 且必须能够接受矢量输入 end else error 输入参数错误 请参考函数帮助 end Ck cotescoeff n Ak b a Ck y Ak fx 2 2 在命令窗口输入 edit cotescoeff 来建立 cotescoeff m 文件 文件中的程序如下 function Ck cotescoeff n for i 1 n 1 k i 1 Ck i 1 n k factorial k factorial n k n quadl t intfun t n k 0 n end 3 在命令窗口输入 edit intfun 来建立 intfun m 文件 文件中的程序如下 function f intfun t n k f 1 for i 0 k 1 k 1 n f f t i end funfun 积分表达式 这里有两种选择 积分表达式 这里有两种选择 1 1 积分函数句柄 必须能够接受矢量输入 比如积分函数句柄 必须能够接受矢量输入 比如 fun x 1 xfun x 1 x 2 x y 2 x y 坐标的离散点 坐标的离散点 第一列为第一列为 x x 第二列为第二列为 y y 必须等距 必须等距 且节点的个数小于且节点的个数小于 9 9 比如 比如 fun 1 8 1 1 8 fun 1 8 1 1 8 如果如果 funfun 的表采用第二种方式 那么只需要输入第一个参数即可 否则还要输入的表采用第二种方式 那么只需要输入第一个参数即可 否则还要输入 a b na