基于MATLAB的双重积分的数值求解(附完整代码与例题)

目录

一. Gauss求积公式

1.1 数学理论

1.2 例题与MATLAB代码

例题1

二. 双重积分问题的数值解

2.1 数学理论

2.2 MATLAB代码与例题

例题2 

例题3 

例题4

例题5

例题6

三. 小结


一. Gauss求积公式

1.1 数学理论

我们知道数值的方法求解积分往往精度会差一点。这个时候Gauss求积公式可以让数值求解积分有较高的代数精度。举例对[-1,1]区间内的Gauss积分:

\int_{-1}^1f(x)dx\approx\sum_{k=0}^nA_kf(x_k)

n取不同值时,Guass求积公式也会不同:

n=1,x_0=x_1=\pm0.5773503

A_0=A_1=1.00000000

n=2,x_0=x_2=\pm0.7745967

A_0=A_2=0.55555556

x_1=0.00000000,A_1=0.88888889

类推到任意的求积区间[a,b],通过变换x=\frac{b-a}{2}t+\frac{b+a}{2}可得如下:

\int_a^bf(x)dx=\frac{b-a}{2}\int_{-1}^1f(\frac{b-a}{2}t+\frac{a+b}{2})dt=\frac{b-a}{2}\sum_{k=0}^nA_kf(\frac{b-a}{2}t+\frac{a+b}{2})

1.2 例题与MATLAB代码

例题1

取n=2,利用Gauss求积公式计算\int_0^1cosxdx

解:

(1)原函数文件

function y=gaussf(x)
y=cos(x);

(2)主运行文件

clc;clear;
gauss2('gaussf',0,1)

function g=gauss2(fun,a,b)
h=(b-a)/2;
c=(a+b)/2;
x=[h*(-0.7745967)+c, c, h*0.7745967+c];
g=h*(0.55555556*(gaussf(x(1))+gaussf(x(3)))+0.88888889*gaussf(x(2)));
end

运行结果:

ans = 0.841471418069758
 

二. 双重积分问题的数值解

2.1 数学理论

矩形区域上的二重积分的数值计算形式:

I=\int_{y_m}^{y_M}\int_{x_m}^{x_M}f(x,y)dxdy

在MATLAB中调用的格式为:

y=dblquad(Fun,x_m,x_M,y_m,y_M)

限定精度的双重积分调用格式为:

y=dblquad(Fun,x_m,x_M,y_m,y_M,\epsilon)

2.2 MATLAB代码与例题

例题2 

求解 

J=\int_{-1}^1\int_{-2}^2e^{\frac{-x^2}{2}}sin(x^2+y)dxdy

解:

MATLAB代码:

clc;clear;
f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');
y=dblquad(f,-2,2,-1,1)

运行结果:

y =1.574493189744943

 在MATLAB中,利用tiled方法也可以计算二重数值积分,调用的函数为quad2d。格式如下:

q=quad2d(fun,a,b,c,d)

此函数可看成近似计算平面区域a\leq x\leq bc(x)\leq y\leq d(x)上,fun(x,y)的积分。fun是函数句柄,c和d可以是标量或者是函数句柄。如果增加计算精度,调用格式如下:

[q,errbnd]=quad2d(fun,a,b,c,d,param1,val1,param2,val2,\ldots)

例题3 

求积分

J=\int_{-\frac{1}{2}}^1\int_{-\sqrt{1-\frac{x^2}{2}}}^{\sqrt{1-\frac{x^2}{2}}}e^{\frac{-x^2}{2}}sin(x^2+y)dydx

解:

MATLAB代码:

clc;clear;

%s数值方法
fh=@(x)sqrt(1-x.^2/2);
fl=@(x)-sqrt(1-x.^2/2); 
f=@(x,y)exp(-x.^2/2).*sin(x.^2+y);
I1=quad2d(f,-1/2,1,fl,fh) % MATLAB自带函数求解一般区域的二重积分

%解析解方法
syms x y  
i2=int(exp(-x^2/2)*sin(x^2+y),y,-sqrt(1-x^2/2),sqrt(1-x^2/2));
I2=int(i2,x,-1/2,1);
I2=vpa(I2)

运行结果:
I1 =0.411929546211955
I2 =0.41192954617629511965175994017601

例题4

计算单位圆域上的积分:

解: 

先把二重积分转化为常见的形式如下:

I=\int_{-1}^1\int^{\sqrt{1-y^2}}_{-\sqrt{1-y^2}} e^{-\frac{x^2}{2}}sin(x^2+y)dxdy

MATLAB代码如下:

clc;clear;

%s数值方法
fl=@(y)-sqrt(1-y.^2); 
fh=@(y)sqrt(1-y.^2); 
f=@(x,y)exp(-x.^2/2).*sin(x.^2+y); 
I1=quad2d(f,-1,1,fl,fh)


%解析解方法
syms x y;  
i2=int(exp(-x^2/2)*sin(x^2+y), x, -sqrt(1-y.^2), sqrt(1-y.^2));
I2=int(i2,y,-1,1);
I2=vpa(I2) 

运行结果:
I1 =0.536860382211891

 
I2 =0.5368603826988078755775938492913
 

最后再介绍一个二重积分进行数值计算的另一个函数Integral2。在MATLAB中调用的格式为如下:

q=integral2(fun,xmin,xmax,ymin,ymax)

如果增加参数限制的话,调用格式如下:

q=integral2(fun,xmin,xmax,ymin,ymax,Name,Value)

例题5

求函数f(x,y)=\frac{1}{(\sqrt{x+y})(1+x+y)}0\leq x\leq10\leq y\leq 1-x三角形区域的积分。

解:

MATLAB代码:

clc;clear;
fun =@(x,y)1./(sqrt(x+y).*(1+x+y).^2);
ymax=@(x)1-x;
q=integral2(fun,0,1,0,ymax)

运行结果:

q =0.285398175390866

例题6

求函数f(x,y)=ax^2+by^2在区域0\leq x\leq5-5\leq y\leq0的积分。其中参数a=3,b=5

解:

MATLAB代码:

clc;clear;
a=3; 
b=5;
fun=@(x,y)a*x.^2+b*y.^2;
q=integral2(fun,0,5,-5,0,'Method','iterated',...
'AbsTol',0,'RelTol',1e-10)

运行结果:q = 1.666666666666667e+03

三. 小结

MATLAB是矩阵和实验室的结合,即矩阵工厂( matrix laboratory),是美国mathworks公司在互联网快速发展的时代下开发的,它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中,将原本极其复杂的问题简单化。

评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

唠嗑!

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值