利用MATLAB求符号微积分

摘要

本文是《科学计算与MATLAB语言》专题七第2小节的学习笔记，如果大家有时间的话，可以去听听课，没有的话，可以看看我的笔记，还是很不错的。本节主要讲了如何利用MATLAB求符号函数的极限、倒数、积分，文中每个代码我都跑过一遍，可以直接复制到MATLAB中运行。

1 符号函数的极限

（1）极限

求符号函数极限的命令为limit，其调用格式为：
limit（f，x，a）
即求函数f关于变量x在a点的极限。

（2）单边极限

limit函数的另一种功能是求单边极限，其调用格式为：
limit（f，x，a，‘right’）
limit（f，x，a，‘left’）
例1 求下列极限。
(1) ${\lim }_{x\to a}\frac{\sqrt[3]{x}-\sqrt[3]{a}}{x-a}$
(2) ${\lim }_{x\to \infty }(1+\frac{1}{n}{)}^n$

syms a m x n ;
f=(x^(1/m)-a^(1/m))/(x-a);
g=(1+1/n)^n;
limit(f,x,a)
limit(g,n,inf)

ans = 
a^(1/m - 1)/m
ans = 
exp(1)

exp(1)即e

2 符号函数的导数

MATLAB中的求导函数为：
diff(f,x,n)
即求函数f关于变量x的n阶导数。参数x的用法同求极限函数limit，可以缺省，默认值与limit相同，n的默认值是1。
例2 求下列函数的导数。
(1) $y=\sqrt{1+e^x}$
(2) $z=\frac{x{e}^y}{y^2},求z_x^{\prime },z_y^{\prime }。$
(1) $y=\sqrt{1+e^x}$

syms x y z;
f=sqrt(1+exp(x));
diff(f)

ans = 
exp(x)/(2*(exp(x) + 1)^(1/2))

(2) $z=\frac{x{e}^y}{y^2},求z_x^{\prime },z_y^{\prime }。$

g=x*exp(y)/y^2;
diff(g,x)
diff(g,y)

ans = 
exp(y)/y^2 
ans = 
(x*exp(y))/y^2 - (2*x*exp(y))/y^3

3 符号函数的积分

（1）不定积分

在MATLAB中，求不定积分的函数是
int()
其常用的调用格式为：
int(f,x)
即求函数f对变量x的不定积分。
例3 求下列不定积分。
(1) $\int (3-x^2{）}^{3}dx$
(2) $\int \frac{5xt}{1+x^2}dt$
(1) $\int (3-x^2{）}^{3}dx$

syms x t;
f=(3-x^2)^3;
int(f)

ans =
- x^7/7 + (9*x^5)/5 - 9*x^3 + 27*x

(2) $\int \frac{5xt}{1+x^2}dt$

g=5*x*t/(1+x^2);
int(g,t)

ans =
(5*t^2*x)/(2*(x^2 + 1))

（2）定积分

在MATLAB中，定积分的计算也使用int（）函数，但调用格式有区别：
int(f,x,a,b)
其中，a、b分别表示定积分的下限和上限。
当函数f关于变量x在闭区间[a,b]可积时，函数返回一个定积分结果。
当a、b中有一个是inf时，函数返回一个广义积分。
当a、b中有一个符号表达式时，函数返回一个符号函数。
例4 求下列定积分。
(1) ${\int }_1^2|1-x|dx$
(2) ${\int }_{-\infty }^{+\infty }\frac{1}{1+x^2}dx$
(3) ${\int }_2^{sinx}\frac{4x}{t}dt$
(1) ${\int }_1^2|1-x|dx$

syms x t;
int(abs(1-x),1,2)

ans =
1/2

(2) ${\int }_{-\infty }^{+\infty }\frac{1}{1+x^2}dx$

int(1/(1+x^2),-inf,inf)

ans = 
pi

(3) ${\int }_2^{sinx}\frac{4x}{t}dt$

int(4*x/t,t,2,sin(x))

ans =
4*x*(log(sin(x)) - log(2))

河道水流量的估算问题
根据实际测量，得到河流某处宽600m，其横截面不同位置某一时刻的水深如下表所示。

x	0	50	100	150	200	250	300	350	400	450	500	550	600
h（x）	4.4	4.5	4.6	4.8	4.9	5.1	5.4	5.2	5.5	5.2	4.9	4.8	4.7

①若此刻水流的流速为0.6m/s，试估计该河流此刻的流量。
②已知x方向[50,60]区间为坡式护岸的下部护脚部分，根据相关堤防设计规范，抛石护岸护脚坡度应缓于1：1.5（正切值），请估计水流冲刷是否已破坏该区域的护脚。

①先拟合出河床曲线，然后进行定积分，计算出河流横截面，即可估计流量。
②根据河床曲线，计算其导函数，并判断相应范围内导函数的取值是否大于1：1.5。

xi=0:50:600; 
yi=[4.4,4.5,4.6,4.8,4.9,5.1,5.4,5.2,5.5,5.2,4.9,4.8,4.7]; 
p=polyfit(xi, yi,3); 
plot(xi, yi,'o', xi, polyval(p, xi)); %同时画出散点图与拟合函数。
syms y x; 
y=poly2sym(p,x); 
s=int(y,x,0,600); 
v=s*0.6; 
eval(v)

ans =
   1.7874e+03

方法一：

xi=0:50:600; 
yi=[4.4,4.5,4.6,4.8,4.9,5.1,5.4,5.2,5.5,5.2,4.9,4.8,4.7]; 
yn=-yi;
p=polyfit(xi, yn,3); %这两句是为了让函数开口朝上，更符合河道真实情况。
plot(xi, yn,'o', xi, polyval(p, xi)); 
syms y x yii; 
y=poly2sym(p,x); 
yii=diff(y,x); 
x=50:60; 
k=eval(yii); 
all(abs(k)<1/1.5)%判断相应范围内导函数的取值是否小于1：1.5。all（i）：若向量i中所有元素非零，结果为1，否则结果为0。

ans =
  logical
   1

方法二：

xi=0:50:600; 
yi=[4.4,4.5,4.6,4.8,4.9,5.1,5.4,5.2,5.5,5.2,4.9,4.8,4.7]; 
yn=-yi;
p=polyfit(xi, yn,3); %这两句是为了让函数开口朝上，更符合河道真实情况。
plot(xi, yn,'o', xi, polyval(p, xi)); %polyval(p, xi)可以求得多项式在xi的函数值并放到p中。
x=50:60; %步长为1
y=polyval(p,x); 
k=diff(y)/1;%diff为差分函数，步长为1，k就是微商，导数。
all(abs(k)<1/1.5) 

ans =
  logical
   1

可以看出，结果为1，也就是说，在[50,60]内,导函数的取值是小于1：1.5，符合要求。

结语

最后欢迎大家点赞👍，收藏⭐，转发🚀，
如有问题、建议，请您在评论区留言💬哦。