B站台湾大学郭彦甫｜MATLAB 学习笔记｜10 数值微积分 Integration-Differentiation

MATLAB学习笔记（10 数值微积分 Integration-Differentiation）

1. 微分 (Differentiation)

1.1 计算多项式：polyval()

example:

画出多项式：

$$9 x^{3}-5 x^{2}+3 x+7$$

​ x 的取值范围为 $-2 \leq x \leq 5$

a = [9,-5,3,7]; %可以看成是原函数的系数矩阵,分别表示x^3,x^2,x^1和x^0对应的系数
x = -2: 0.01: 5;
f = polyval(a,x);% y = polyval(a,x) 计算多项式 a 在 x 的每个点处的值

plot(x,f,'LineWidth', 2);
xlabel('x'); ylabel('f(x)');
set(gca, 'FontSize', 14)

1.2 计算多项式微分：polyder()

example:

原函数为：

$$f\left( x \right) =5x^4-2x^2+1$$

求：

（1）求导数 $f'(x)$

p=[5 0 -2 0 1]; %原函数的系数矩阵
polyder(p)	% 得到一个向量，为微分后的系数矩阵

（2）求 $f'(7)$

polyval(polyder(p),7)	
% polyval(polyder(p),7)	表示求微分后的函数在 7 处的取值。

exercise (P7)

题目：

画出多项式：

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

和该函数的导数函数，x 的取值范围为 $-2\le x\le 1$

解：

a = [5 -7 5 10];
b = [0 4 12 -3];
c = conv(a,b); % conv()为卷积函数，对表示系数的数列进行卷积即为多项式相乘
x = -2: 0.01: 1;
f = polyval(c,x); % 得到原函数 f(x)
s = polyval(polyder(c),x); % 得到求导后的函数 s(x)

hold on
plot(x,f,'k','LineWidth', 2);
plot(x,s,'r','LineWidth', 2);
xlabel('x'); ylabel('f(x)');
legend({'原函数f(x)','导数f''(x)'},'Location','northwest');
title("exercise");
hold off

2. 积分 (Integration)

2.1 计算多项式积分：ployint()

给出以下函数：

$$f\left( x \right) =5x^4-2x^2+1$$

计算：

（1）$\int{f\left( x \right) dx}+3$

a = [5 0 0 -1 0 1];
f = polyint(a,3); % 计算多项式 a 的积分，其中常数项为 3 
% q = polyint(p,k) 使用积分常量 k 返回 p 中系数所表示的多项式积分

（2）$\int{f\left( 7 \right) dx}+3$

F = polyval(polyint(a,3),7);    % 计算积分在 x=7 处的取值

3. 数值微分 (Numerical Differerntiation)

（用于非多项式函数，比如正弦函数 $sin(x)$）

数学基础：

在连续的情况下，导数的定义为：

$$f^{\prime}\left(x_{0}\right)=\lim _{h \rightarrow 0} \frac{f\left(x_{0}+h\right)-f\left(x_{0}\right)}{h}$$

其中，$h$ 表示 $x$ 非常小的变化。

3.1 计算微分用：diff()

计算离散点的导数：

$$f'\left( x \right) =\frac{\Delta \text{y}}{\Delta \text{x}}=\frac{y_1-y_2}{x_1-x_2}$$

例子：

计算点 (1,5) 和点 (2,7) 构成的线之间的斜率：

x = [1 2]; y = [5 7];
slope = diff(y)./diff(x);
% diff() 计算系数矩阵中元素间前后的差值

问题：

给出原函数 $f(x)=sin(x)$，找出 $f'(x)$ 在 $x_0=\pi/2$ 的导数，且 $h=0.1$.

解：

利用导数定义：

$$f^{\prime}\left(x_{0}\right)=\lim _{h \rightarrow 0} \frac{f\left(x_{0}+h\right)-f\left(x_{0}\right)}{h}$$

计算：

$$f'\left( \frac{\pi}{2} \right) =\lim_{h\rightarrow 0}\frac{f\left( \frac{\pi}{2}+h \right) -f\left( \frac{\pi}{2} \right)}{h}$$

代码如下：（取步进 h=0.1）

x0 = pi/2; h = 0.1;
x = [x0 x0+h];
y = [sin(x0) sin(x0+h)];
m = diff(y)./diff(x)

exercise (P13)

题目：

原函数为 $f(x)=sin(x)$，写出一个程序：根据变量 $h$ 的不同，找到 $f'(x_0)$ 在 $x_0=\pi/2$ 的误差。

（理想的 $f'(\frac{\pi}{2})=0$）

解：

x0 = pi/2;
for n=0:6
    h = 0.1/10^n;
    x = [x0 x0+h];
    y = [sin(x0) sin(x0+h)];
    f(n+1) = diff(y)./diff(x);
end

h	error of $f'(x_0)$
0.1	-0.049958347219743
0.01	-0.004999958333474
0.001	-4.999999583255583e-04
0.0001	-4.999999969613196e-05
0.00001	-5.000000413669099e-06
0.000001	-5.000444503323075e-07
0.0000001	-4.996003607896202e-08

例子：

画出 $f(x)=sin(x)$ 及其导数在 $x=[0,2\pi]$ 范围内的图像：

clear 
clf
h = 0.1;
x = 0:h:2*pi;
y = sin(x);
m = diff(y)./diff(x);

hold on
plot(x,y,'k','LineWidth',2);
plot(x(1:end-1),m,'r','LineWidth',2);   
% 计算微分得到的系数向量的维度（m的维度）=原函数的系数向量的维度（x的维度）-1
title('exercise');
legend('sin(x)','sin''(x)','Location','northeast');
hold off

exercise (P17)

题目：

给出原函数 $f\left( x \right) =e^{-x}\sin \left( x^2/2 \right)$，画出近似导数 $f'$ 在 $h=0.1,\ 0.01, \ 0.001$ 处的曲线。

解：

clear
clf
for n = 0:1:2
    
    h = 0.1/10^n;
    x = 0:h:2*pi;
    hold on
    if n == 0
        y0 = exp(-x) .* sin(x.^2./2);
        m0 = diff(y0)./diff(x); % m为导数
        plot(x(1:end-1),m0,'b','LineWidth',1);
        
    elseif n == 1
        y1 = exp(-x) .* sin(x.^2./2);
        m1 = diff(y1)./diff(x); % m为导数
        plot(x(1:end-1),m1,'g','LineWidth',1);
        
    elseif n == 2
        y2 = exp(-x) .* sin(x.^2./2);
        m2 = diff(y2)./diff(x); % m为导数
        plot(x(1:end-1),m2,'r','LineWidth',1);
        
    end
    hold off
end

title('exercise');
legend('$f''\left( x \right)$ (h=0.1)',...
    '$f''\left( x \right)$ (h=0.01)',...
    '$f''\left( x \right)$ (h=0.001)',...
        'Location','northeast','Interpreter','latex');

（放大后的曲线如下，三条线都已经画出来了）

3.2 二阶和三阶导数

给出原函数为 $f(x) = x^3$ ，求 $f'(x)$ 和 $f''(x)$ 在 $-2\le x \le 2$ 的取值。

x = -2:0.005:2; y = x.^3;
m = diff(y)./diff(x);
m2 = diff(m)./diff(x(1:end-1));	% 每做一次微分，系数矩阵中的元素数量-1

4. 数值积分 (Numerical Integration)

数学基础：

计算定积分：

$$s=\int_{a}^{b} f(x) d(x) \approx \sum_{i=0}^{n} f\left(x_{i}\right) \int_{a}^{b} L_{i}(x) d x$$

使用求面积的方法：（使用一个有限的点集去近似整数）

4.1 基本求面积方法

4.1.1 中点法（Midpoint rule（零阶近似））

Example：

计算

$$A=\int\limits_0^2{4x^3dx=}x^4\mid _{0}^{2}=\left( 2 \right) ^4-\left( 0 \right) ^4=16$$

h = 0.05; x = 0:h:2;
midpoint = (x(1:end-1)+x(2:end))./2;
y = 4*midpoint.^3;
s = sum(h*y)	% 计算矩形面积用长 h 乘以高 h

>> class10

s =

   15.9950

4.1.2 梯形方法 （Trapezoid rule（一阶近似））

%方法一
h = 0.05; x = 0:h:2; y = 4*x.^3;
s = h*trapz(y)	% 不明白 trapz() 的用法

>> class10

s =

   16.0100

%方法二
h = 0.05; x = 0:h:2; y = 4*x.^3;
trapezoid = (y(1:end-1)+y(2:end))/2;	% 梯形公式
s = h*sum(trapezoid)

>> Class14

s =

   16.0100

4.2 二阶方法：$\frac{1}{3}$ 辛普森公式

example：

计算：

$$A=16\int\limits_0^2{4x^3dx=}x^4\mid _{0}^{2}=\left( 2 \right) ^4-\left( 0 \right) ^4=16$$

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))

4.2.1 前面三种方法对比

4.3 计算积分用：integral()

integral() 采用全局自适应计算面积方法和默认的误差容忍度对函数进行数值积分

example：

计算：

$$A=\int\limits_0^2{\frac{1}{x^3-2x^2-5}dx}$$

y = @(x) 1./(x.^3-2*x-5);	% 使用指针@ 表示函数，不能直接使用表达式
integral(y,0,2)