【台大郭彦甫】Matlab入门教程超详细学习笔记七：数值微积分（附PPT链接）

---

前言

【台大郭彦甫】PPT链接：https://pan.baidu.com/s/1VXdy3HCBPexMK0csyehiZA 提取码：2io1
matlab官方帮助文档：https://ww2.mathworks.cn/help/

---

微分

• 函数f(x)的导数写成f‘(x) 或者 df(x)，表示函数f(x)相对于x的变化率。在几何上，f‘(xo)表示点xo与曲线切线方向的变化量，也就是斜率。

一、多项式微积分

1. 多项式计算

polyval()

• 多项式微分表达式如下：
  

• matlab如何表示多项式？使用行向量

• y = polyval(p,x)计算多项式 p 在 x 的每个点处的值。参数 p 是长度为 n+1 的向量，其元素是 n 次多项式的系数（降幂排序）。

eg.

a = [9,-5,3,7];
x = -2:0.01:5;
f = polyval(a,x);
plot(x,f,'LineWidth', 2);
xlabel('x');
ylabel('f(x)');
set(gca,'FontSize', 14)

2. 多项式微分

polyder()

• k = polyder(p) 返回 p 中的系数表示的多项式的导数
• k = polyder(a,b) 返回多项式 a 和 b 的乘积的导数
• [q,d] = polyder(a,b) 返回多项式 a 和 b 的商的导数

eg.

• 对多项式求微分

p=[5 0 -2 0 1];
polyder(p)

结果：

• 计算x=7处的微分值

polyval(polyder(p),7)

练习

• 提示：conv()卷积和多项式乘法
• w = conv(u,v) 返回向量 u 和 v 的卷积。如果 u 和 v 是多项式系数的向量，对其卷积与将这两个多项式相乘等效。

x=-2:0.01:1;
a1=[5,-7,5,10]
a2=[4,12,-3];
a=conv(a1,a2);      %计算两多项式相乘所得多项式系数
y=polyval(a,x);
a_=polyder(a);      %计算f(x)微分式的系数
y_=polyval(a_,x);
plot(x,y,'--b',x,y_,'r','linewidth',2);
legend('f(x)','f''(x)');

注意：

• plot()中设置线条宽度的'linewidth'参数，对前面所画两条线都起作用；
• 添加图例时，在字符串中显示单引号，打两个单引号'f''(x)'即显示一个单引号。

3. 多项式积分

polyint()

• 多项式积分表达式如下：
  
• q = polyint(p,k) 使用积分常量 k 返回 p 中系数所表示的多项式积分。
• q = polyint(p) 假定积分常量 k = 0。

eg.

• 对多项式求微分，指定常数项为3

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

结果：

• 计算x=7处的积分值

polyval(polyint(p, 3),7)

结果：

二、数值微积分

1. 数值微分法

diff()

• 数值微分表达式如下：
  
  
  Y = diff(X) 计算相邻元素之间的差分

eg.

x = [1 2 5 2 1];
diff(x)

结果：

x = [1 2];y = [5 7];
slope = diff(y)./diff(x)

结果：

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

结果：

• h越小误差越小

练习

x0 = pi/2;
h = 0.1;
for i = 1:1:7
    x = [x0 x0+h];
    y = [sin(x0) sin(x0+h)];
    error = diff(y)./diff(x);
    A = ['h=',num2str(h),' error=',num2str(error)];
    disp(A)
    h = h.*0.1;
end

h	error of f’(x)
0.1	-0.049958
0.01	-0.0050
0.001	-0.0005
0.0001	-5.0000e-05
0.00001	-5.0000e-06
0.000001	-5.0000e-07
0.0000001	-4.9960e-08

• 如何找到0~2$\pi$区间上的$f^{\prime }$
  策略：
  ①在间隔[0，2𝜋]中创建数组
  ②步骤是ℎ
  ③计算这些点的𝑓′
  

h = 0.5;
x = 0:h:2*pi;
y = sin(x);
y_wf = diff(y)./diff(x);
plot(x,y);
x(length(x)) = [];
hold on
plot(x,y_wf);
legend('sin(x)',"sin'(x)");
set(gca, 'FontSize', 15)

• g = colormap(lines);其实是把一个256x3的颜色矩阵[RGB]赋值给g，然后方便给各曲线设置不同颜色
• g(i,:)即为256x3的颜色矩阵中第i行对应的颜色

g = colormap(lines);
hold on;
for i=1:4
    x = 0:power(10, -i):pi;
    y = sin(x);
    y_ds = diff(y)./diff(x);
    plot(x(1:end-1), y_ds, 'Color',  g(i,:));
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;

练习

g = colormap(lines);
hold on;
for i=1:3
    x = 0:power(10, -i):2*pi;
    y = exp(-x).*sin(x.^2./2);
    y_ds = diff(y)./diff(x);
    plot(x(1:end-1), y_ds, 'Color',  g(i,:),'LineWidth',2);
end
hold off;
set(gca, 'xlim', [0, 2*pi]);
set(gca, 'ylim', [-0.3, 0.3]);
set(gca, 'FontSize', 18);
set(gca, 'XTick', 0:pi/2:2*pi);
set(gca, 'XTickLabel', {'0',  '\pi/2','\pi','3\pi/2', '2\pi'});
h =  legend('h=0.1','h=0.01','h=0.001');
set(h,'FontName', 'Times New  Roman');
box on;

2. 高阶微分法

• 一次微分：m=diff(f(x))./diff(x)
• 二次微分：m2=diff(m)./diff(x)
• 二阶导数𝑓′′和三阶导数可以用类似的方法得到

eg.$$f(x)=x^3$$

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);
xlabel('x', 'FontSize', 18);
ylabel('y', 'FontSize', 18);
legend('f(x) = x^3','f''(x)','f''''(x)');
set(gca, 'FontSize', 18);

3. 数值积分法

1. 中点法则（零级近似）—— sum()
   
   

eg.

思路：

• $x=[x_1{x}_2{x}_3{x}_4\dots x_{end}]$
• $x(1:end-1)=[x_1{x}_2{x}_3{x}_4\dots x_{end-1}]$
  $x(2:end)=[x_2{x}_3{x}_4\dots x_{end}]$
• $\frac{x(1:end-1)+x(2:end)}{2}=[\frac{x_1+x_2}{2}\frac{x_2+x_3}{2}\frac{x_3+x_4}{2}...\frac{x_{end-1}+x_{end}}{2}]$

eg.

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，提高精度

eg.

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

结果：

2. 梯形法则（一阶近似）—— trapz()
   
   
   eg.
   

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

结果：

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

结果：

3. 二阶法则：1/3辛普森法则【比较精准】

eg.

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. 三种方法的比较：

三、回顾Function Handles(@)

Function Handles (@)

• Function Handles即函数句柄，是一种表示函数的 MATLAB数据类型。
• 函数句柄的典型用法是将函数传递给另一个函数。

eg.

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

• 先将上述代码保存为：xy_plot.m文件，然后运行下面代码

xy_plot(@sin,0:0.01:2*pi);

xy_plot(@cos,0:0.01:2*pi);

xy_plot(@exp,0:0.01:2*pi);

四、直接计算积分和微分

1. 数值积分：integral()

eg.

y = @(x) 1./(x.^3-2*x-5);
integral(y,0,2)

2. 二重积分：integral2()

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

3. 三重积分integral3()

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

---

总结

以上就是第七节的内容，本部分介绍了matlab的数值微积分部分。