matlab一维矩形积分,玩转matlab之一维 gauss 数值积分公式及matlab源代码
在数值分析中,尤其是有限元刚度矩阵、质量矩阵等的计算中,必然要求如下定积分: $$ I=\int_a^b f(x)dx $$学好gauss积分也是学好有限元的重要基础,学过高等数学的都知道,手动积分能把人搞死(微笑脸),而且有些函数还不存在原函数,使用原始的手动算出原函数几乎是不现实的。因此非常有必要学习数值积分,简单讲就是近似计算,只要这个近似值精确度高和稳定性好就行。Gauss积分公式就是这么一个非常好用的工具。本文介绍高斯积分公式的使用以及简单的数值算例。
标准区间
先考虑特殊情况,对于一般区间呢?待会会处理这个问题。 $$ I=\int_{-1}^1 f(x)dx $$ 不加证明的直接给出gauss公式如下:详情参阅任何一本数值分析书都有详细的证明过程: $$ I=\int_{-1}^1 f(x)dx=\Sigma_{i=1}^n A_if(x_i) $$ 其中$A_i$称作权,$x_i$称作 gauss 点。 下面的问题就是如何选择$n,A_i,x_i$。
理论表明n个点的Gauss公式代数精度为$2n-1$,其选择如下表,(这里仅仅举1-4个点情况,实际使用的时候一般2点或者3点的精度已经完全够了)更多积分点可参考 gauss表.
gauss点个数 $n$
gauss 点 $x_i$
权重 $A_i$
精度
1
$x_1$=0
$A_1$=2
1
2
$x_{1,2}=\pm1/\sqrt{3}$
$A_1=A_2=1$
3
3
$x_1=-\sqrt{3/5}$
$x_2=0$
$x_3=\sqrt{3/5}$
$A_1=5/9$
$A_2=8/9$
$A_3=5/9$
5
4
$x_{1}=-\sqrt{\dfrac{15+2\sqrt{30}}{35}}$
$x_{2}=-\sqrt{\dfrac{15-2\sqrt{30}}{35}}$
$x_{3}=\sqrt{\dfrac{15-2\sqrt{30}}{35}}$
$x_{4}=\sqrt{\dfrac{15+2\sqrt{30}}{35}}$
$A_1=\frac{90-5\sqrt{30}}{180}$
$A_2=\frac{90+5\sqrt{30}}{180}$
$A_3=\frac{90+5\sqrt{30}}{180}$
$A_4=\frac{90-5\sqrt{30}}{180}$
7
一般区间
$$ I=\int_a^b f(x)dx $$
根据上面的讨论情况,可知只要做变换(相当于换元积分一样) $$ 令\quad x=\frac{b+a+(b-a)s}{2},\ 则\quad dx = \frac{b-a}{2}ds. $$ 那么有$s\in[-1,1]$,于是即可使用标准区间公式如下: $$ I = \int_a^bf(x)dx=\int_{-1}^1f(\frac{b+a+(b-a)s}{2})\times\frac{b-a}{2}ds\ =\frac{b-a}{2}\Sigma_{i=1}^mA_if(\frac{b+a+(b-a)s_i}{2}) $$ 上述公式中的$A_i$即为表格中的权重,$s_i$即为上表中对应的gauss点,代入公式即可计算积分值。
数值实验
[toc] 所有实验在MATLAB2018a版本下完成。(建议安装新版本,因为很多函数在新版中已经优化了或者改名字了,比如老版本积分函数quad 新版已经改为integral,只不过目前quad函数还是可以使用的,将来会被删除)。
我们取2个函数做实验,分别计算出其gauss积分值再与matlab自带的函数 integral 计算结果作比较,实验模型是: $$ 计算 \quad I= \int_1^2 f(x)dx $$
实验一
取函数 $$ f(x)=lnx, \quad 即自然对数函数以e为底. $$ 使用matlab函数 integral 计算得到: $I= 0.386294361119891$。 使用gauss积分的matlab计算结果为:
高斯点数 m
积分值 $I_m$
误差norm($I_m-I$)
2
0.386594944116741
3.01E-04
3
0.386300421584011
6.06E-06
4
0.386294496938714
1.36E-07
5
0.386294364348948
3.23E-09
实验二
取函数 $$ f(x)=\dfrac{x^2+2x+1}{1+(1+x)^4}, $$ 使用matlab函数 integral 计算得到: $I= 0.161442165779443$。 使用gauss积分的matlab计算结果为:
高斯点数 m
积分值 $I_m$
误差norm($I_m-I$)
2
0.161394581386268
4.76E-05
3
0.161442818737102
6.53E-07
4
0.161442196720137
3.09E-08
5
0.161442166345131
5.66E-10
总结
随着gauss点m的个数增多,精度在逐渐提高,但是要注意的是,gauss点取得多的话,计算量也会增大,只是因为我们计算的问题规模比较小,所以感觉不到而已。
另外可以看到2点3点的gauss公式的精度已经很高了,说明并不需要取太多的点,而在实际计算中,比如有限元的计算中,也仅仅取2点或者3点gauss积分就完全足够。
下节预告
下次介绍gauss积分的二维公式使用以及matlab数值实验,欢迎有问题与我交流,偏微分方程,矩阵计算,数值分析等问题,我的qq 群 315241287
matlab代码
clc;clear;
% using 2 3 4 5 points compute the integral
% x \in [a,b]
%
%% test
a=1; b = 2;
fun = @(x) log(x);
% fun = @(x) 2*x./(1+x.^4);
% fun = @(x) exp(-x.^2/2);
% fun = @(x) (x.^2+2*x+1)./(1+(1+x).^4);
%% setup the gauss data
for gauss = 2:5
if gauss == 2
s=[-1 1]/sqrt(3);
wt=[1 1];
fprintf('*************************** 2 points gauss *******')
elseif gauss == 3
s = [-sqrt(3/5) 0 sqrt(3/5)];
wt = [5 8 5]/9;
fprintf('*************************** 3 points gauss *******')
elseif gauss == 4
fprintf('*************************** 4 points gauss *******')
s = [ -sqrt((15+2*sqrt(30))/35), -sqrt((15-2*sqrt(30))/35), ...
sqrt((15-2*sqrt(30))/35), sqrt((15+2*sqrt(30))/35)];
wt = [ (90-5*sqrt(30))/180, (90+5*sqrt(30))/180,...
(90+5*sqrt(30))/180, (90-5*sqrt(30))/180];
elseif gauss == 5
fprintf('*************************** 5 points gauss *******')
s(1)=.906179845938664 ; s(2)=.538469310105683;
s(3)=.0; s(4)=-s(2) ; s(5)=-s(1);
wt(1)=.236926885056189 ; wt(2)=.478628670499366;
wt(3)=.568888888888889 ; wt(4)=wt(2) ; wt(5)=wt(1);
end
%%
% 区间变换到 s \in[-1,1]
s = (b-a)/2*s+(b+a)/2;
jac = (b-a)/2;% dx = jac * ds
f = fun(s);
f = wt.* f .* jac;
format long
exact = integral(fun,a,b);
comp = sum(f)
fprintf('the error is norm(comp-exact)=%10.6e\n\n\n',norm(comp-exact))
end
fprintf('\n\n********* matlab built-in function ''integral''*********\n')
exact = integral(fun,a,b)
format short
matlab一维矩形积分,玩转matlab之一维 gauss 数值积分公式及matlab源代码相关推荐
- 玩转matlab之一维 gauss 数值积分公式及matlab源代码
目录 标准区间 一般区间 数值实验 实验一 实验二 总结 下节预告 matlab代码 在数值分析中,尤其是有限元刚度矩阵.质量矩阵等的计算中,必然要求如下定积分: \[ I=\int_a^b f(x) ...
- matlab的积分公式,玩转matlab之一维 gauss 数值积分公式及matlab源代码
释放双眼,带上耳机,听听看~! 目录 在数值分析中,尤其是有限元刚度矩阵.质量矩阵等的计算中,必然要求如下定积分: \\[ I=\\int_a^b f(x)dx \\]学好gauss积分也是学好有限元 ...
- 玩转 matlab 之一维 gauss 数值积分公式及matlab源代码
文章目录 标准区间 一般区间 数值实验 实验一 实验二 总结 下节预告 matlab代码 在数值分析中,尤其是 有限元刚度矩阵.质量矩阵等的计算中,必然要求如下定积分: I=∫abf(x)dxI=\i ...
- 运用数学软件matlab求无穷积分,matlab积分的计算及其简单应用论文.doc
积分的计算及其简单应用 摘要:本文简要的概述了MATLAB 在高等数学中积分的计算及应用:利用MATLAB 中符号积分和数值积分的命令,计算定积分和不定积分.同时,也可以通过这些命令来解决一些实际问题 ...
- 基于MATLAB的特殊函数积分
目录 前言 (一)振荡函数的积分 (二)反常(广义)积分 1. 无界函数的反常积分 2. 无穷区间上的反常积分 一. quadgk()函数在MATLAB中的运用 二. 基于MATLAB的特殊函数积分 ...
- matlab 积分函数曲线,Matlab之函数积分 | 学步园
Matlab之函数积分 一元函数(一重)积分: 求一元函数积分有quad函数,quadl函数,int函数: quad和quadl: quad和quadl两个函数,他们使用不同的 例:求 的积分: f= ...
- matlab 微分命令 求导,Matlab微分和积分
第六讲 Matlab 微分和积分 理论介绍:微分.有限差分.积分.离散求和 软件求解:函数及常见注意事项 一.一元函数导数与微分 Matlab 由命令函数diff 来完成求导运算,调用格式为:diff ...
- 普朗克公式matlab,用MATLAB实现普朗克函数积分的快捷计算.pdf
用MATLAB实现普朗克函数积分的快捷计算 维普资讯 ≮ 曩≯ : 一 l2 量 ii!iiiiiiliii! 文章编号: 1672-8785(2008)04-0012-03 用MATLAB实现普朗克 ...
- MATLAB | 全网唯一,使用MATLAB绘制矩形树状图
绘制效果 全网唯一这四个字我都快说腻了,请叫我绘图小天才,又双叒叕写了一个工具函数发到了MATHWORKS,矩形树状图主要用于直观展示各个元素的分类和占比. 编写不易点个赞叭~~ 基本使用 需要准备一 ...
最新文章
- 将服务器置于最终用户附近可解决性能问题?—Vecloud微云
- pfile文件怎么恢复格式_回收站清空的文件怎么恢复?值得收藏的恢复方法
- Android 长按照片保存 工具类
- linux删除、读取文件原理
- devexpress gridcontrol 内置导航栏,双击后才修改数据
- linux include 编译,linux-如何使用OpenSSL include编译.c文件?
- 四面体的表面积_边长为正四面体的表面积是()、;、;、;、。
- [转载] Java异常:选择Checked Exception还是Unchecked Exception?
- pyspark指定schema
- 年轻人“躺平”的生活方式,引起不少争议
- 精灵3P+Pix4D简单航测详细应用教程
- python合并pdf_一个用于合并pdf的简单Python脚本
- power supply surges detected
- c语言 步进电机 程序,步进电机加速启动C语言程序
- php设为首页代码,JavaScript
- 互联网行业的裁员潮;程序员到35岁是个坎儿!
- 不一样的“中国速度”,数据可视化交通运输大屏,带你见证中国高铁
- android反编译去壳,安卓apk查壳工具,逆向反编译必备
- 软考高级-系统架构师-案例分析-数据库真题考点汇总
- CSS3弹性盒子布局
热门文章
- 阿里云商标注册申请进度查询方法
- stm32 步进电机控制,S曲线加减速,匀加速运动
- 在带有双硬盘的Windows10系统上安装Ubuntu16.04系统
- ~杂记(3):los_dispatch.s和startup.s的作用分析
- 国外计算机领域的ei期刊,2008 EI 收录国外英文期刊(计算机类)
- RPO和RTO是什么?
- 蘑菇街Java后台面试总结
- 计算时间复杂度--(简单版)
- docker:error pulling image configuration
- 美国计算机专业四年毕业率达多少,卡内基梅隆大学毕业率知多少