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

piday123

# m4 z" E! f/ r5 m目录

• 总述
• 函数说明
• 应用举例
• 函数实现
  & ?! k* M# ~; g" @" @5 {

" X: _7 D5 F8 U9 _4 w- ~9 N& ^! [
总述
1 i' r6 f( ~' C5 ?# {2 {' V  @1 z3 T$ s

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

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

0 a/ D' Y% s9 L; `$ w$ F" F函数说明

% R% X2 n) m, i5 A4 k, t5 i+ e

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

2 \! {/ o* `) M# \7 m& Y应用举例

问题： 求函数

的1~4阶导数, 并验证误差。

0 J6 B2 _! ]5 a! a0 h4 r

代码如下：

  o" `  Q: C. \

• % // 输入函数，并求解析解，并代入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)))
  

, Z- k" y  K! I& y* n
$ {4 a: K% d, L2 Y- U) \. l' Z; p4 U

不同阶的导数图像如下：

$ L: V+ R" `4 ]8 Z1 j

定量地分析误差时， 考虑到计算得出的4阶导数向量， 其长度比原始对照向量f4短， 所以两个向量取同样多点进行比较， 就可以得出数值方法的相对误差最大值为

, 亦即0.035%。 由此可见， 这里的数值方法还是很精确的。

! t3 v4 W0 Q3 a9 n* k函数实现
' g. X/ j- t  I% \" P- f: L

• 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;
  

0 b! X3 n# a0 L7 C$ O2 o

xiaogegepcb

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