数值微积分

前言一、多项式微积分1. 多项式计算2. 多项式微分3. 多项式积分二、数值微积分1. 数值微分法2. 高阶微分法3. 数值积分法三、回顾Function Handles(@)四、直接计算积分和微分1. 数值积分：integral()2. 二重积分：integral2()3. 三重积分integral3()总结

---

前言

【台大郭彦甫】PPT链接：/s/1VXdy3HCBPexMK0csyehiZA 提取码：2io1

matlab官方帮助文档：/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:7x = [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

如何找到0~2π\piπ区间上的f′f'f′

策略：

①在间隔[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 onplot(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:4x = 0:power(10, -i):pi;y = sin(x);y_ds = diff(y)./diff(x);plot(x(1:end-1), y_ds, 'Color', g(i,:));endhold 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:3x = 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);endhold 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)=x3f(x)=x^3f(x)=x3

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. 数值积分法

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

eg.

思路：

x=[x1x2x3x4…xend]x = [ x_1   x_2   x_3   x_4   …   x_{end} ]x=[x1​ x2​ x3​ x4​ … xend​]x(1:end−1)=[x1x2x3x4…xend−1]x(1:end−1) = [ x_1   x_2   x_3   x_4   …   x_{end-1} ]x(1:end−1)=[x1​ x2​ x3​ x4​ … xend−1​]

x(2:end)=[x2x3x4…xend]x(2:end) = [x_2   x_3   x_4   …   x_{end}]x(2:end)=[x2​ x3​ x4​ … xend​]x(1:end−1)+x(2:end)2=[x1+x22x2+x32x3+x42...xend−1+xend2]\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}]2x(1:end−1)+x(2:end)​=[2x1​+x2​​2x2​+x3​​2x3​+x4​​...2xend−1​+xend​​]

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)

结果：

梯形法则（一阶近似）—— 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 xy = 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的数值微积分部分。