基于MATLAB的数值积分问题求解（附完整代码与例题）

唠嗑！

已于 2023-12-30 16:23:37 修改

阅读量1.8w

点赞数 25

分类专栏： MATLAB 文章标签：线性代数 matlab 图论数学建模计算机网络

于 2022-03-22 19:38:10 首次发布

本文链接：https://blog.csdn.net/forest_LL/article/details/123624830

版权

MATLAB 专栏收录该内容

45 篇文章 106 订阅

订阅专栏

一. 由给定数据进行梯形求积分

1.1 数学理论

梯形求积的图形理解可看如下图：

把函数分解成无数个这种类似的梯形，如下：

积分形式：

$I=\int_a^b{f(x)}dx=\sum_{i=1}^N\int_{x_i}^{x_{i+1}}f(x)dx=\sum_{i=1}^N\Delta f_i$

梯形面积形式：

$S=\frac{1}{2}[\sum_{i=1}^{N-1}(y_{i+1}+y_i)(x_{i+1}-x_i)]=\frac{1}{2}\lbrace\sum_{i=1}^{N-1}[(y_{i+1}-y_i)+2y_i](x_{i+1}-x_i)\rbrace$

上式中向量:

$x={[x_1,x_2,\ldots,x_N]}^T,y={[y_1,y_2,\dots,y_N]}^T$

梯形面积对应的MATLAB代码：

sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2

MATLAB自带有函数：trapz() 来计算数值梯形积分。调用格式为：

trapz(x,y)

来看几道例题。

1.2 例题与MATLAB代码

例题一

用梯形法求 $x\in(0,\pi)$ 区间内，函数 $sin(x),cos(x),sin(\frac{x}{2})$ 的定积分值。

解：

MATLAB代码：

clc;clear;
x1=[0:pi/30:pi]';
y=[sin(x1) cos(x1) sin(x1/2)];
x=[x1 x1 x1];


S1=sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2 %直接求解

S2=trapz(x1,y) %调用函数求解

运行结果：

S1 =1.998171961343654 0.000000000000000 1.999543052990808

S2 =1.998171961343654 0.000000000000000 1.999543052990808

结果分析：两种方法得出的结果完全一致

例题二

用定步长方法求解积分 $\int_0^{\frac{3\pi}{2}}cos(15x)dx$

解：

MATLAB代码如下：

clc;clear;

%画图
x=[0:0.01:3*pi/2]; %这样赋值能确保3*pi/2点包含在内
y=cos(15*x);
plot(x,y)

%求取理论值
syms x;
A=int(cos(15*x),0,3*pi/2)

%数值方法
h0=[0.1,0.01,0.001,0.0001,0.00001,0.000001]; %随着步距h的减小，计算精度逐渐增加
v=[];
for h=h0
    x=[0:h:3*pi/2];
    y=cos(15*x);
    I=trapz(x,y);
    v=[v;h,I,1/15-I]; %1/15是积分的理论值
end
format long;
v

运行结果：

A =1/15

结果的解释：第一列为不同的步长，第二列为对应求出的积分值，第三列为与标准积分结果的差值。

二. 单变量数值积分问题求积

2.1 数学理论

首先引入辛普森(Simpson)公式。Simpson方法求解 $[x_i,x_{i+1}]$ 上的积分如下：

$\Delta f_i\approx \frac{h_i}{12}[f(x_i)+4f(x_i+\frac{h_i}{4})+2f(x_i+\frac{h_i}{2})+4f(x_i+\frac{3h_i}{4})+f(x_i+h_i)]$

上式子中 $h_i=x_{i+1}-x_i$

形成图形解释，如下：

MATLAB中可以调用函数直接计算如下：

y=quad(Fun,a,b)  y=quadl(Fun,a,b) % 求定积分
y=quad(Fun,a,b,c) y=quadl(Fun,a,b,c)
%两种函数均为变步长
%Fun为函数的字符串变量
%限定精度的定积分求解，默认精度为10－6。
%quadl函数使用的算法是自适应Lobatto算法其精度和速度均远高于quad函数

2.2 例题与MATLAB代码

例题三

求解 $erf(1.5)=\frac{2}{\sqrt{\pi}}\int_0^{1.5} e^{-t^2}dt$

代码如下：

clc;clear;

%运用符号工具箱
syms x;
y0=vpa(int(2/sqrt(pi)*exp(-x^2),0,1.5),60)

format long %16位精度
f=@(x)2/sqrt(pi)*exp(-x.^2);
y=quadl(f,0,1.5,1e-20) %设置高精度

error=abs(y-y0)

运行结果：

y0 =0.966105146475310713936933729949905794996224943257461473285748

y =0.966105146475311

error =0.000000000000000064025228489138917324084325807202

例题四

求解I并画出f(x)。

$I=\int_0^4f(x)dx$ ，其中f(x)定义为如下:

解：

（1）画图

代码如下：

x=[0:0.01:2, 2+eps:0.01:4,4];
y=exp(x.^2).*(x<=2)+80./(4-sin(16*pi*x)).*(x>2);
y(end)=0;
x=[eps, x];
y=[0,y]; %为减少视觉上的误差，对端点与间断点（有跳跃）进行处理。

fill(x,y,'g')

运行结果：

（2）求积分

代码如下：

clc;clear;
f=@(x)exp(x.^2).*(x<=2)+80./(4-sin(16*pi*x)).*(x>2);

%直接调用quadl求解
I1=quadl(f,0,4)

%解析解
syms x;
I2=vpa(int(exp(x^2),0,2)+int(80/(4-sin(16*pi*x)),2,4))

运行结果：

I1 =57.764450169467679
I2 =57.76445012505301033331523538518

以下介绍一个新的MATLAB可调用的函数：integral数值积分

q=integral(fun,xmin,xmax,Name,Value)

例题五

计算函数 $f(x)=e^{-x^2}{(lnx)}^2$ 的积分， $x\in[0,Inf)$ 。（此题为广义积分）

解：

MATLAB代码如下：

clc;clear;
fun=@(x)exp(-x.^2).*log(x).^2;
q=integral(fun,0,Inf)

运行结果：q =1.947522220295560

例题六

计算参数化函数 $f(x)=\frac{1}{x^3-2x-c}$ 的积分， $x\in[0,2],c=5$ 。

MATLAB代码如下：

clc;clear;
c=5;
fun=@(x,c)1./(x.^3-2*x-c);
q=integral(@(x)fun(x,c),0,2)

运行结果：q = -0.460501533846733

再补充举一个例子：

代码：

clc;clear;
fun=@(x)log(x);
format long
q1=integral(fun,0,1)
q2=integral(fun,0,1,'RelTol',0,'AbsTol',1e-12)

运行结果：

q1 =-1.000000010959678

q2 =-1.000000000000010

唠嗑！

关注

25
点赞
踩
167

收藏

觉得还不错? 一键收藏
打赏
1
评论
基于MATLAB的数值积分问题求解（附完整代码与例题）

第一列为不同的步长，第二列为对应求出的积分值，第三列为与标准积分结果的差值。以下介绍一个新的MATLAB可调用的函数：integral数值积分。首先引入辛普森(Simpson)公式。Simpson方法求解。运行结果：q = -0.460501533846733。运行结果：q =1.947522220295560。来计算数值梯形积分。：两种方法得出的结果完全一致。求解I并画出f(x)。用定步长方法求解积分。
复制链接

扫一扫