玩转matlab之一维 gauss 数值积分公式及matlab源代码

目录

  • 标准区间
  • 一般区间
  • 数值实验
    • 实验一
    • 实验二
  • 总结
  • 下节预告
  • 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点,代入公式即可计算积分值。

数值实验

所有实验在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

总结

  1. 随着gauss点m的个数增多,精度在逐渐提高,但是要注意的是,gauss点取得多的话,计算量也会增大,只是因为我们计算的问题规模比较小,所以感觉不到而已。
  2. 另外可以看到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

原文地址:https://www.cnblogs.com/sunzhenwei/p/10847059.html

时间: 2024-11-09 03:08:04

玩转matlab之一维 gauss 数值积分公式及matlab源代码的相关文章

MATLAB/Excel-如何将Excel数据导入MATLAB中

在使用MATLAB对矩阵进行数据处理时,为了方便编辑与修改,常常需要先将数据录入到Excel中,然后再将其导入到MATLAB中参与矩阵运算.本文在MATLAB 2013a和Office 2013环境下向大家演示如何将Excel数据导入到MATLAB中,其他版本的MATLAB.OFFICE方法大同小异,一起来看一下 工具/原料   Excel数据文件(格式xls或xlsx) MATLAB 7.x + 方法/步骤     将待导入的矩阵结构的数据录入Excel中,录入时注意行列要跟原矩阵一一对应  

【Matlab编程】生日快乐歌(显示歌词)—matlab版

clear all A4=440;%标准音A4 不同的曲调音调不同scale的取值范围不同 pt=44100;p0=pt/2;%频率 scale=A4/2^(9/12)*2.^((-12:11)/12);%这里可以调节音调高低,eg:改变式子中的-12:11为0:23 map=[1 3 5 6 8 10 12 13 15 17 18 20 22 24 25];%音符,这个需要看曲谱编码 score=[5 5 6 5 8 7 5 5 6 5 9 8 5 5 12 10 8 7 6 6 11 11

Matlab与C/C++联合编程之Matlab以MEX方式调用C/C++代码(三)

最近写了个Matlab程序,好慢呐……所以开始学习Matlab与C/C++混合编程.下面写了个测试代码,显示一个Double类型矩阵中的元素. 源代码 #include "mex.h" void displaySubscript( const mxArray *pArray, mwSize index ); // 入口函数void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) {   

Matlab高级教程_第二篇:Matlab相见恨晚的模块_02_并行运算-2

1 MATLAB并行计算-从个人桌面到远程集群和云(陈伟/魏奋)视频摘录笔记 https://cn.mathworks.com/videos/parallel-computing-with-matlab-92865.html 2 数据.硬件和算法,MATLAB发展的方向 3 MATLAB并行计算工具的介绍 内嵌多线程(隐式) --MATLAB 内核函数和图像处理工具箱 --矩阵操作(linear algebra,fft,filter,etc) --无须代码修改 并行计算产品(显式) --Para

一维信号频谱图仿真——matlab

程序1: %在MATLAB中是用连续信号在等时间间隔点的样值来近似地表示连续信号的,当采样时间间隔足够小时,这些离散的采样值就能较好地近似出连续信号,matlab中连续信号的显示实际上还是离散信号的显示,只是取样点特别 %多的时候,用线连接起来,显示出来的图形就比较圆滑,接近连续信号:如果取样点特别少,连接起来就会变成折线: clear all; %这个其实可以没有,只不过以前出过问题,现在就加上! N=1024; %这个是你举得信号的点数,随便你了 fs=16000; %这个是抽样频率,记得要

Matlab的XTickLabel中数值带下标

%axis为'x'或'y',分别表示更改x或y刻度 %ticks是字符cell function settick(axis,ticks) n=length(ticks); tkx=get(gca,'XTick');tky=get(gca,'YTick'); switch axis case 'x' w=linspace(tkx(1),tkx(end),n); set(gca, 'XTick', w, 'XTickLabel', []);%刷新刻度,去掉刻度值 yh=(14*w(1)-w(end)

Matlab高级教程_第二篇:Matlab相见恨晚的模块_02_并行运算-利用GPU并行执行MATLAB程序

1 MATLAB原文: 如果所有你想使用的函数支持GPU,你能够使用gpuArray把输入的数据传输到GPU,也能够唤起gather命令把传输值GPU的数据回收. 2 通过gpuDevice命令观察当前电脑的GPU设备 >> gpuDevice ans = CUDADevice (具有属性): Name: 'GeForce GT 430' % GPU设备的型号 Index: 1 % 当前GPU设备的编号 ComputeCapability: '2.1' % 计算能力 SupportsDoubl

Matlab与C/C++联合编程之Matlab以MEX方式调用C/C++代码(二)

如果我有一个用C语言写的函数,实现了一个功能,如一个简单的函数: double add(double x, double y) { return x + y; } 现在我想要在Matlab中使用它,比如输入: >> a = add(1.1, 2.2) 3.3000 要得出以上的结果,那应该怎样做呢? 解决方法之一是要通过使用MEX文件,MEX文件使得调用C函数和调用Matlab的内置函数一样方便.MEX文件是由原C代码加上MEX文件专用的接口函数后编译而成的.可以这样理解,MEX文件实现了一种

matlab学习笔记第四章——统计和MATLAB编程介绍

1.柱状图: >> x = [55,63,69,70,75,78,82,84,85,88,90,96,100]; >> y = [1,2,1,6,4,7,2,1,3,2,4,2,1]; >> bar(x,y) 2.我们可以使用barh命令产生水平的柱状图. 3.通过调用mean函数,MATLAB会告诉我们一组数据的均值是多少. 4.