matlab实现数值微分（diff_ctr函数）

清疚 2022-11-25 00:59 255阅读 0赞

### 目录 ###

*   *   *  总述
         *  函数说明
         *  应用举例
         *  函数实现

### 总述 ###

如果**已知**函数表达式，可以通过`diff(）`函数求取各阶导数解析解的方法，并得出结论，高达100阶的导数也可以用MATLAB语言在几秒钟的时间内直接求出。

如果函数表达式**未知**，只有实验数据，在实际应用中经常也有求导的要求，这样的问题就不能用前面的方法获得问题的解析解。要求解这样的问题，需要引入数值算法得出所需问题的解。由于在MATLAB语言中没有现成的数值微分函数，所以本文将介绍一种数值微分算法——**中心差分方法**。

### 函数说明 ###

function [dy,dx] = diff_ctr(y,Dt,n)
    %diff_ctr
    %中心差分算法实现数值微分
    %  调用格式：
    %    [d_y, d_x] = diff_ctr(y,Dt,n)
    %  其中，y为给定的等间距的实测数据构成的向量， Dt为自变量的间距，n为所需的导数阶次。 
    %  向量d_y为得出的导数向量， 而d_x为相应的自变量向量。注意这两个向量的长度比y短。
    %
    % Examples:
    %  求函数y=sin(x)/(x^2+4*x+3)的1~4阶导数
    % MATLAB求解语句：
    %  h=0.05; x=0:h:pi; syms x1; 
    %  f=sin(x1)/(x1^2+4*x1+3); y=subs(f,x1,x);
    %  [y1,dx1]=diff_ctr(y,h,1); subplot(221), plot(dx1,y1); 
    %  [y2,dx2]=diff_ctr(y,h,2); subplot(222), plot(dx2,y2); 
    %  [y3,dx3]=diff_ctr(y,h,3); subplot(223), plot(dx3,y3); 
    %  [y4,dx4]=diff_ctr(y,h,4); subplot(224), plot(dx4,y4);
    
    % 与解析解对比验证：
    % syms x1; 
    % f=sin(x1)/(x1^2+4*x1+3); 
    % yy1=diff(f);   f1=subs(yy1,x1,x); 
    % yy2=diff(yy1); f2=subs(yy2,x1,x); 
    % yy3=diff(yy2); f3=subs(yy3,x1,x); 
    % yy4=diff(yy3); f4=subs(yy4,x1,x);
    % % 求四阶导数向量的范数（相对误差）：
    % norm(double((y4-f4(4:60))./f4(4:60)))

### 应用举例 ###

**问题：** 求函数  y = s i n x x 2 + 4 x + 3 y=\\frac\{sin x\}\{x^2+4x+3\} y=x2\+4x\+3sinx 的1~4阶导数, 并验证误差。

**代码如下：**

% // 输入函数，并求解析解，并代入x向量得出精确解。
    h=0.05; x=0:h:pi; syms x1; 
    f=sin(x1)/(x1^2+4*x1+3); 
    yy1=diff(f); f1=subs(yy1,x1,x); 
    yy2=diff(yy1); f2=subs(yy2,x1,x); 
    yy3=diff(yy2); f3=subs(yy3,x1,x); 
    yy4=diff(yy3); f4=subs(yy4,x1,x);
    %// 比较不同阶的导数
    y=subs(f,x1,x); 
    [y1,dx1]=diff_ctr(y,h,1); subplot(221), plot(x,f1,dx1,y1,':'); 
    [y2,dx2]=diff_ctr(y,h,2); subplot(222), plot(x,f2,dx2,y2,':'); 
    [y3,dx3]=diff_ctr(y,h,3); subplot(223), plot(x,f3,dx3,y3,':'); 
    [y4,dx4]=diff_ctr(y,h,4); subplot(224), plot(x,f4,dx4,y4,':')
    %// 定量分析误差
    norm(double((y4-f4(4:60))./f4(4:60)))

不同阶的导数**图像**如下：  
![函数图像][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70]  
定量地分析误差时， 考虑到计算得出的4阶导数向量， 其长度比原始对照向量`f4`短， 所以两个向量取同样多点进行比较， 就可以得出数值方法的相对误差最大值为 3.5 × 1 0 − 4 3.5 \\times 10^\{-4\} 3.5×10−4, 亦即  0.035 % 0.03 5\\% 0.035% 。 由此可见， 这里的数值方法还是很精确的。

### 函数实现 ###

function [dy,dx] = diff_ctr(y,Dt,n)
    y1=[y 0 0 0 0 0 0];
    y2=[0 y 0 0 0 0 0];
    y3=[0 0 y 0 0 0 0];
    y4=[0 0 0 y 0 0 0];
    y5=[0 0 0 0 y 0 0];
    y6=[0 0 0 0 0 y 0];
    y7=[0 0 0 0 0 0 y];
    switch n
        case 1
            dy = (-y1+8*y2-8*y4+y5)/12/Dt;
        case 2
            dy = (-y1+16*y2-30*y3+16*y4-y5)/12/Dt^2;
        case 3
            dy = (-y1+8*y2-13*y3+13*y5-8*y6+y7)/8/Dt^3;
        case 4
            dy = (-y1+12*y2-39*y3+56*y4-39*y5+12*y6-y7)/6/Dt^4;
    end
    dy = dy(5+2*(n>2):end-4-2*(n>2));
    dx = ([2:length(dy)+1]+(n>2))*Dt;

--------------------

此函数源文件可前往下面网址下载：

[diff\_ctr.m下载通道][diff_ctr.m]

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70]: /images/20221124/f7eea995fff045bd806f31b19dc82d8d.png
[diff_ctr.m]: https://download.csdn.net/download/weixin_43964993/12693884