玩转 matlab 之二维 gauss 数值积分公式使用及 matlab 源代码(1)-常量区间
目录
继续上一篇一维gauss积分的讨论,本文讨论二维gauss数值积分公式的使用,并给出数值实验。
一. 标准区间
按照一维情况类似方式,首先给出如下二重积分公式:
\[
I = \int_{-1}^1dx\int_{-1}^1 f(x,y)dy=\Sigma_{k=1}^m\Sigma_{l=1}^mA_{kl}f(x_k,y_l).
\]
这是标准的二维gauss积分 \(m\times m\) (即 x 轴和 y 轴分别取 m 个点对应\(m^2\)个点)公式的表达式,其中,\(x_{kl}\) 称作gauss积分点,\(A_{kl}\)称作\(x_{kl}\) 点处的权重。\(m=1,2,3\) 时候草图如下所示:
下面给出部分gauss积分点和权重对应列表
gauss点个数 \(m^2(m\times m)\) | gauss 点 \((x_i,y_i)\) | 权重 \(A_i\) | 精度(2m-1) |
---|---|---|---|
1 (\(1\times 1\)) | \(x_1=y_1=0\) | \(A_1\)=4 | 1 |
4 (\(2 \times 2\)) | \(x_{1}=-1/\sqrt{3},y_{1}=-1/\sqrt{3}\) \(x_{2}=1/\sqrt{3},y_{2}=-1/\sqrt{3}\) \(x_{3}=1/\sqrt{3},y_{3}=1/\sqrt{3}\) \(x_{4}=-1/\sqrt{3},y_{4}=1/\sqrt{3}\) |
\(A_1=A_2=A_3=A_4=1\) | 3 |
9 (\(3 \times 3\)) 令\(gpt=\sqrt{3/5}\) |
\(x_1=-gpt,y_1=-gpt\) \(x_2=gpt,y_2=-gpt\) \(x_3=gpt,y_3=gpt\) \(x_4=-gpt,y_4=gpt\) \(x_5=0,y_5=-gpt\) \(x_6=gpt,y_6=0\) \(x_7=0,y_7=gpt\) \(x_8=-gpt,y_8=0\) \(x_9=0,y_9=0\) |
\(A_1=A_2=A_3=A_4=25/81\) \(A_5=A_6=A_7=A_8=40/81\) \(A_9=64/81\) |
5 |
一般二重积分近似值也就是使用 \(2\times 2,\quad 3\times 3\) 公式就完全足够了。所以不需要太多的点,徒增计算量。因此就可以得到gauss积分的坐标表达式为:
\[
I = \int_{-1}^1dx\int_{-1}^1 f(x,y)dy=\Sigma_{i=1}^{m^2}A_{i}f(x_i,y_i).
\]
二. 一般区间
对于一般区间,先考虑区间端点为常量情况(下一节介绍区间为变量情况),即
\[
I = \int_{a}^bdx\int_{c}^d f(x,y)dy.
\]
其中 \(a,b,c,d\) 都是已知常数。与一维情况类似,只需要做变量变换,于是\([s,t]\in[-1,1]\times [-1,1]\)(通俗讲就是换元法)
\[
令 \quad x =\frac{b+a+(b-a)s}{2},\qquad 则\quad dx=\frac{b-a}{2}ds,\\
\quad y =\frac{d+c+(d-c)t}{2},\qquad 则\quad dy=\frac{d-c}{2}dt.\\
记\qquad jac = \frac{(b-a)(d-c)}{4},\qquad 则 \quad dxdy=jac\times dsdt
\]
于是二重积分就变成了
\[
I = \int_{a}^bdx\int_{c}^d f(x,y)dy
=jac\times\int_{-1}^1ds\int_{-1}^1f(\frac{b+a+(b-a)s}{2},\frac{d+c+(d-c)t}{2})dt.\\=jac\times\Sigma_{i=1}^{m^2}A_{i}f(\frac{b+a+(b-a)s_i}{2},\frac{d+c+(d-c)t_i}{2}).
\]
其中 \((s_i,t_i)\) 即为上表中的 gauss 节点,对应的权重因子为 \(A_i\)。
三. 数值实验(没有编程的计算是不完整的)
使用matlab2018a 计算结果,并且与matlab自带函数 integral2 计算的结果进行比较给出误差。
算例如下:
\[
计算定积分 I = \int_{a}^bdx\int_{c}^d f(x,y)dy
\]
其中,\(a=1.4,b=2,c=1,d=1.5,f(x,y)=ln(x+2*y), ln\)是以e为底对数函数。使用matlab的integral2 函数计算结果为\(I =0.429554527548275\).
自己编程计算结果如下:
高斯点数 \(m^2(m\times m)\) | 积分值 \(I_m\) | 误差norm(\(I_m-I\)) |
---|---|---|
4(\(2\times 2\)) | 0.429556088022242 | 1.56E-06 |
9(\(3\times 3\)) | 0.429554531152490 | 3.60E-09 |
四. 总结和下节预告
- 从实验数据可以发现,二重gauss数值积分使用\(2\times 2\) 4个点的情况下,结果已经很准确了,达到了1e-6的误差,所以在实际计算中,一般使用4点或者9点的gauss公式就可以满足要求了。
- 下一节介绍当积分区间在变系数下的二重gauss公式的计算方法
- 欢迎与我交流,数值分析,矩阵计算,PDE数值解等,QQ群 315241287
五. matlab源代码
clc;clear; % compute int_a^b [int_c)^d f(x,y)] % (x,y) \in [a,b] X [c,d] %% setup the integral interval and gauss point and weight a = 1.4; b = 2; c = 1; d =1.5; fun=@(x,y) log(x+2*y); fprintf('***********************************************\n') for gauss = 2:3 % m points rule in 2 dimensional case if gauss == 2 fprintf('******* 2X2 points gauss rule result *******') gpt=1/sqrt(3); s(1) = -gpt; t(1) = -gpt; s(2) = gpt; t(2) = -gpt; s(3) = gpt; t(3) = gpt; s(4) = -gpt; t(4) = gpt; wt = [1 1 1 1]; elseif gauss == 3 gpt=sqrt(0.6); fprintf('******* 3X3 points gauss rule *******') s(1) = -gpt; t(1) = -gpt; wt(1)=25/81; s(2) = gpt; t(2) = -gpt; wt(2)=25/81; s(3) = gpt; t(3) = gpt; wt(3)=25/81; s(4) = -gpt; t(4) = gpt; wt(4)=25/81; s(5) = 0.0; t(5) = -gpt; wt(5)=40/81; s(6) = gpt; t(6) = 0.0; wt(6)=40/81; s(7) = 0.0; t(7) = gpt; wt(7)=40/81; s(8) = -gpt; t(8) = 0.0; wt(8)=40/81; s(9) = 0.0; t(9) = 0.0; wt(9)=64/81; end %% 区间变换到 [-1,1] X [-1,1] jac = (b-a)*(d-c)/4; x = (b+a+(b-a)*s)/2; y = (d+c+(d-c)*t)/2; f = fun(x,y); comp = wt(:) .* f(:) .* jac;%无论一个向量是行还是列,写成x(:)都会变成列向量 format long comp = sum(comp) exact = integral2(fun,a,b,c,d); fprintf('the error is norm(comp-exact)=%10.6e\n\n',norm(comp-exact)) end fprintf('******************************************\n') fprintf('matlab built-in function ''integral2''\n') exact format short
- 玩转matlab之一维 gauss 数值积分公式及matlab源代码
- 数值积分之Gauss求积法五点公式
- 使用MATLABD数值法计算定积分或反常积分
- 使用RTL-SDR和Matlab Simulink玩转软件无线电(二)
- (数值计算方法)matlab的使用
- MATLAB学习笔记(八)——MATLAB数值积分与微分
- 一个使用MATLAB手动求二维曲线交点的例子
- 应朋友们要求,Hello China的源代码下载资源积分取消,欢迎大家下载使用
- 【数学建模】MATLAB数值积分与微分
- 相机标定简介与MatLab相机标定工具箱的使用(未涉及原理公式推导)
- (数值计算方法)matlab的使用
- MATLAB学习笔记:数值积分
- 使用RTL-SDR和Matlab Simulink玩转软件无线电(二十一)
- 应朋友们要求,Hello China的源代码下载资源积分取消,欢迎大家下载使用
- 使用RTL-SDR和Matlab Simulink玩转软件无线电(二十四)
- 韩顺平_PHP程序员玩转算法公开课(第一季)09_使用栈完成高级计算器(1)_学习笔记_源代码图解_PPT文档整理
- 使用RTL-SDR和Matlab Simulink玩转软件无线电(十一)
- 使用RTL-SDR和Matlab Simulink玩转软件无线电(八)
- 使用RTL-SDR和Matlab Simulink玩转软件无线电(一)
- 数值分析基础工具使用Matlab绘制双曲线