微积分计算

多项式的微积分

多项式微分

matlab中使用向量来表示多项式.即在向量中只存储多项式对应次数项的系数即可. 比如多项式f(x)=x3−2x−5f(x)=x^{3}-2 x-5在MATLAB中可以表示为 p=[1 0 -2 -5 ];

使用polyval函数可以获得多项式函数在某些点处的值, 得到这些函数值后就可以通过plot函数来绘制多项式函数的图像.比如想要绘制9x3−5x2+3x+79 x^{3}-5 x^{2}+3 x+7在−2≤x≤5-2 \leq x \leq 5处的函数图像, 可以通过执行下面的代码实现.

1
2
3
4
5
a = [9,-5,3,7]; x = -2:0.01:5;
f = polyval(a,x); % f中存储了多项式a在点x处的数值
plot(x,f,'LineWidth', 2);
xlabel('x'); ylabel('f(x)');
set(gca,'FontSize', 14)

画出的图形如下:

image-20200714130747862

多项式函数的求导可以使用polyder函数来实现, 既然多项式表示成向量形式, 那其导函数自然也是向量形式.

对于求函数f(x)=5x4−2x2+1f(x)=5 x^{4}-2 x^{2}+1 在x=7处的导数的问题, 可以通过polyval函数和polyder函数共同实现.

1
2
3
p=[5 0 -2 0 1];
polyder(p)
polyval(polyder(p),7)

其结果为:

image-20200714131330480

可以验证其结果是正确的的.

练习: 绘制多项式函数图像及其导数图像

f(x)=(5x3−7x2+5x+10)(4x2+12x−3)f(x)=\left(5 x^{3}-7 x^{2}+5 x+10\right)\left(4 x^{2}+12 x-3\right)

其范围在−2≤x≤1-2 \leq x \leq 1

提示: 可能用到的函数conv()

1
2
3
4
5
6
7
8
9
10
11
p=[5 -7 5 10];
q=[4 12 -3];
f=conv(p,q);
df=polyder(f);
x=linspace(-2,1);
hold on
plot(x,polyval(f,x),'r','LineWidth', 2);
plot(x,polyval(df,x),'--b','LineWidth', 2);
hold off
set(gca,'FontSize', 14);
legend('f(x)','f^{\prime}(x)');

参考效果图:

image-20200714133602909

多项式积分

多项式函数的积分使用polyint来实现, 需要注意的是, 对于一个多项式函数, 其不定积分结果有无穷多个, 因此在使用polyint函数来求积分时, 不仅要给定一个被积函数, 还需要给出积分结果的常数项的值.

下述代码实现对多项式f(x)=5x4−2x2+1f(x)=5x^4-2x^2+1的积分, 令常数项为3, 并计算积分后函数在x=7处的函数值.

1
2
3
p=[5 0 -2 0 1];
polyint(p, 3)
polyval(polyint(p, 3),7)

计算结果:

image-20200714134126870

数值微积分

对于一般的曲线函数, 比如f(x)=sinxf(x)=sinx等, matlab计算其导数的方法使用的导数定义法, 即通过导数定义式计算每个点的导数值.

f′(x0)=lim⁡h→0f(x0+h)−f(x0)hf^{\prime}\left(x_{0}\right)=\lim _{h \rightarrow 0} \frac{f\left(x_{0}+h\right)-f\left(x_{0}\right)}{h}

数值微分

对于这样一个复杂的导数定义式, 我们可以分步对它进行计算. 先计算出Δy\Delta y和Δx\Delta x, 在计算其比值.

matlab中提供diff函数可以用来计算向量中相邻两个数的差值. 例如:

1
2
3
4
x = [1 2 5 2 1];
diff(x);
% 结果为:
ans = [1 3 -3 1]

不难发现, 使用diff函数求得的结果值, 就是我们需要的一系列的Δx\Delta x的值. 因此, 想要求曲线的导数值, 只需用diff(y)除以diff(x)即可.

1
2
3
4
5
% 计算f(x)=sin(x)在x等于pi/2处的导数值 定义式中取h=0.1
x0 = pi/2; h = 0.1;
x = [x0 x0+h];
y = [sin(x0) sin(x0+h)];
m = diff(y)./diff(x) % 注意应使用 ./

结果为:

image-20200714140035328

由于, sin(x)在x=π/2x=\pi/2处的导数值是0, 不是-0.05, 因此, 这种计算导数的方式是一种逼近式的计算, 我们可以通过调整h的大小, 来不断的提高精度.

练习:

计算f(x)=sin(x)f(x)=sin(x)在x=π/2x=\pi/2处的导数值, 分别取h= 0.1,0.01,0.001,0.0001,0.00001,0.000001,0.0000001观察结果.

参考代码:

1
2
3
4
5
6
7
8
9
h=[0.1 0.01 0.001 0.0001 0.00001 0.000001 0.0000001];
x0=pi/2;
m=zeros(size(h));
for ii=1:size(h,2)
	x=[x0 x0+h(ii)];
	y=[sin(x0) sin(x0+h(ii))];
	m(ii)=diff(y)./diff(x);
end
m

执行结果:

image-20200714141726667

可以看到, 误差逐渐变小, 由于精度提高, 后三位数显示的都是0; 可以通过执行format long来得到更加精确的数值.

image-20200714141933322

相同的思想可以计算出f(x)=sin(x)f(x)=sin(x)在区间[0,2π][0 ,2\pi]上的导函数.

1
2
3
4
h = 0.5;
x = 0:h:2*pi;
y = sin(x);
m = diff(y)./diff(x);

进而可以绘制出导函数的图像.我们可以尝试使用不同的步长 h 来绘制导函数的图像.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
g = colormap(lines);
hold on;
for i=1:4
	x = 0:power(10, -i):pi; % power()函数是乘方函数
	y = sin(x);
	m = diff(y)./diff(x);
	plot(x(1:end-1), m, 'Color', g(i,:)); % 注意绘制图像时调整x的向量长度
end
hold off;
set(gca,'XLim', [0, pi/2]);
set(gca,'YLim', [0, 1.2]);
set(gca,'FontSize', 18);
set(gca,'FontName','symbol');
set(gca,'XTick', 0:pi/4:pi/2);
set(gca,'XTickLabel', {'0', 'pi/4', 'pi/2'});
h = legend('h=0.1','h=0.01','h=0.001','h=0.0001');
set(h,'FontName', 'Times New Roman'); box on;

结果图形:

image-20200714145003158

这看起来是两条线, 是因为随h的减小, 精度越来越高导致有一定的重合, 可以通过点击放大发现存在着不同的曲线.

练习: 绘制函数f(x)=e−xsin⁡(x2/2)f(x)=e^{-x} \sin \left(x^{2} / 2\right)的导数图像, 分别采用h=0.1,0.01,0.001.

1
2
3
4
5
6
7
8
9
10
11
g = colormap(lines);
hold on;
for i=1:3
	x = 0:power(10, -i):2*pi; % power()函数是乘方函数
	y = exp(-x).*sin(x.*x/2);
	m = diff(y)./diff(x);
	plot(x(1:end-1), m, 'Color', g(i,:)); % 注意绘制图像时调整x的向量长度
end
hold off;
set(gca,'XTick', 0:pi/4:2*pi);
legend('h=0.1','h=0.01','h=0.001');

参考结果图:

image-20200714151211905

根据高阶导数的定义, 高阶导数只不过是其低一阶导函数的导数. 利用这一思想, 可以计算函数的二阶导数, 三阶导数等等.

1
2
3
4
5
6
7
8
x = -2:0.005:2; y = x.^3;
m = diff(y)./diff(x);
m2 = diff(m)./diff(x(1:end-1));
plot(x,y,x(1:end-1),m,x(1:end-2),m2); % 求二阶导数时, x的向量长度应该再减少一个.
xlabel('x','FontSize', 18);
ylabel('y','FontSize', 18);
legend('f(x) =x^3','f^{\prime}(x)','f^{\prime \prime}(x)');
set(gca,'FontSize', 18);

image-20200714151912166

数值积分

根据定积分的几何意义可知, 求一个函数的定积分, 就是求函数图像与坐标轴和积分边界所围成的图形的面积. matlab中使用的是定积分的定义去进行计算, 即"大化小, 常代变"的思想. 把需要计算分面积分割成若干小部分, 计算这些部分图形的面积和.

image-20200714152449291

在图形分割的过程中, 可以使用不同的分割策略, 即使用不同的可计算的多边形面积代替曲边梯形的面积.这里介绍中点规则, 梯形规则和辛普森规则

中点规则

其原理是使用分割边界的中点处的函数值作为分割小矩形的高.即

∫x0x3f(x)dx≈hf0+hf1+hf2=h∑i=0n−1fi\int_{x_{0}}^{x_{3}} f(x) d x \approx h f_{0}+h f_{1}+h f_{2}=h \sum_{i=0}^{n-1} f_{i}

原理图如下:

image-20200714153131470

根据此原理, 计算

A=∫024x3dx=x4∣02=(2)4−(0)4=16A=\int_{0}^{2} 4 x^{3} d x=\left.x^{4}\right|_{0} ^{2}=(2)^{4}-(0)^{4}=16 

1
2
3
4
5
h = 0.05;
x = 0:h:2;
midpoint = (x(1:end-1)+x(2:end))./2;
y = 4*midpoint.^3;
s = sum(h*y)

上述代码的计算结果为:

image-20200714153528442

结果与正确结果16之间存在一定误差, 可以通过减小步长h的方式 来更加细分, 以提高计算精度.

梯形规则

梯形规则的原理是 不再使用小矩形代替曲边梯形, 而是使用直边梯形代替曲边梯形. 即,

∫x0x3f(x)dx≈hf0+f12+hf1+f22+hf2+f32\int_{x_{0}}^{x_{3}} f(x) d x \approx h \frac{f_{0}+f_{1}}{2}+h \frac{f_{1}+f_{2}}{2}+h \frac{f_{2}+f_{3}}{2}

原理示意图:

image-20200714153819868

使用此方法再次计算

A=∫024x3dx=x4∣02=(2)4−(0)4=16A=\int_{0}^{2} 4 x^{3} d x=\left.x^{4}\right|_{0} ^{2}=(2)^{4}-(0)^{4}=16 

1
2
3
h = 0.05; x = 0:h:2; y = 4*x.^3;
trapezoid = (y(1:end-1)+y(2:end))/2;
s = h*sum(trapezoid)

可以直接使用trapz函数来计算trapezoid = (y(1:end-1)+y(2:end))/2;上述代码等价为

1
2
h = 0.05; x = 0:h:2; y = 4*x.^3;
s = h*trapz(y)

执行结果均为:

image-20200714154552434

辛普森规则

辛普森规则的原理推导比较复杂, 这里给出其原理的结论:

∫x0x2f(x)dx≈h3(f0+4f1+f2)\int_{x_{0}}^{x_{2}} f(x) d x \approx \frac{h}{3}\left(f_{0}+4 f_{1}+f_{2}\right)

这个方法一次计算两个步长, 也就是说第一次计算x0x_0到x2x_2, 下一个参与计算的单元为 x2x_2到x4x_4,原理图如下:

image-20200714154951051

使用此规则计算这个积分

A=∫024x3dx=x4∣02=(2)4−(0)4=16A=\int_{0}^{2} 4 x^{3} d x=\left.x^{4}\right|_{0} ^{2}=(2)^{4}-(0)^{4}=16 

1
2
3
h = 0.05; x = 0:h:2; y = 4*x.^3;
s = h/3*(y(1)+2*sum(y(3:2:end-2))+...
4*sum(y(2:2:end))+y(end)) % 开头结尾只加一次, 中间偶数序号每个计算两次, 奇数序号每个计算四次

计算结果为:

image-20200714160122819

这种方法是这三种方法中计算精度最高的方法.

下图中比较了这三种计算方式的区别:

image-20200714160344702

可以理解为, 矩形规则是在用0次曲线逼近, 梯形规则是在用1次曲线逼近, 而辛普森规则使用二次曲线逼近.

内置积分函数

首先复习一下matlab的函数句柄, 以及使用方式.

函数句柄类似于C++中的指针, 它是指向了一个函数. 以下面代码为例介绍函数句柄.

1
2
3
4
5
6
function [y] = xy_plot(input,x)
% xy_plot receives the handle of a function
% and plots that function of x
y = input(x); plot(x,y,'r--');
xlabel('x'); ylabel('function(x)');
end

使用此函数绘图应该传入一个函数句柄, 而不是函数的名称, 正确的使用方式为:

1
2
3
xy_plot(@sin,0:0.01:2*pi);
xy_plot(@cos,0:0.01:2*pi);
xy_plot(@exp,0:0.01:2*pi);

matlab中计算积分的函数是integral函数, 参数为被积函数和上下限. 此处传入函数应该传入其函数句柄.

1
integral(@sin,0,pi)

结果是:

image-20200714161808030

也可以传入我们自定义函数.下面给出计算这个积分的代码.

A=∫024x3dx=x4∣02=(2)4−(0)4=16A=\int_{0}^{2} 4 x^{3} d x=\left.x^{4}\right|_{0} ^{2}=(2)^{4}-(0)^{4}=16 

1
2
f=@(x) 4*x.^3;
integral(f,0,2)

结果为:

image-20200714162101518

计算二重积分和三重积分使用函数integral2和integral3, 下面给出使用示例.

计算: f(x,y)=∫0π∫π2π(y⋅sin⁡(x)+x⋅cos⁡(y))dxdyf(x, y)=\int_{0}^{\pi} \int_{\pi}^{2 \pi}(y \cdot \sin (x)+x \cdot \cos (y)) d x d y

1
2
f = @(x,y) y.*sin(x)+x.*cos(y);
integral2(f,pi,2*pi,0,pi)

计算: f(x,y)=∫−11∫01∫0π(y⋅sin⁡(x)+z⋅cos⁡(y))dxdydzf(x, y)=\int_{-1}^{1} \int_{0}^{1} \int_{0}^{\pi}(y \cdot \sin (x)+z \cdot \cos (y)) d x d y d z

1
2
f = @(x,y,z) y.*sin(x)+z.*cos(y);
integral3(f,0,pi,0,1,-1,1)

资料链接

参考视频 参考讲义