基于MATLAB的数值积分问题求解

您所在的位置：网站首页 › matlab求积分的函数function › 基于MATLAB的数值积分问题求解

基于MATLAB的数值积分问题求解

2023-12-16 07:19| 来源: 网络整理| 查看: 265

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

例题一

例题二

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

例题三

例题四

例题五

例题六

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

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

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

积分形式： $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)

来看几道例题。

例题一

用梯形法求 $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

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

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

首先引入辛普森(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函数例题三

求解 $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).*(x2); y(end)=0; x=[eps, x]; y=[0,y]; %为减少视觉上的误差，对端点与间断点（有跳跃）进行处理。 fill(x,y,'g')

运行结果：

（2）求积分

代码如下：

clc;clear; f=@(x)exp(x.^2).*(x2); %直接调用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

【本文地址】

基于MATLAB的数值积分问题求解

基于MATLAB的数值积分问题求解

今日新闻

推荐新闻