找回密码
 注册
关于网站域名变更的通知
查看: 701|回复: 1
打印 上一主题 下一主题

matlab实现数值微分(diff_ctr函数)

[复制链接]

该用户从未签到

跳转到指定楼层
1#
发表于 2021-1-28 10:21 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式

EDA365欢迎您登录!

您需要 登录 才可以下载或查看,没有帐号?注册

x
本帖最后由 piday123 于 2021-1-28 10:25 编辑
3 Z! M3 H, d2 g$ {1 G
" B4 _4 f! K$ n0 y0 r$ C, E% u& `目录
  • 总述
  • 函数说明
  • 应用举例
  • 函数实现
    4 f& I8 A& E! Z
$ u* X0 Q* U2 v8 r

/ a+ X$ }7 K+ e( Q/ Q8 I总述
. \; A. Z6 e" T) N) X. o6 t8 U3 @1 [: v

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

9 d. s/ o' l$ I/ `- L5 o% S" I

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


0 r& v! o8 X, n# t9 |% i9 i- Z函数说明, K& W4 b: ?( }5 e2 r8 X" O, p
. Y1 }; f) }' G
  • 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)))
      w8 M; }& h* A" S5 T
8 c1 Z* r. A! z6 B
/ X/ \$ K% ?1 n6 Y2 u
: U6 g; m5 J. k8 O& E3 \. D& L8 g
应用举例

问题: 求函数

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


' x% P. S2 T/ ?2 [( u  _: Z

代码如下:


% d2 x7 u' _/ u
  • % // 输入函数,并求解析解,并代入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)))
    3 y8 u7 t  Z+ S" s/ z

% C( l$ T" b. x- b& s5 V6 @/ S3 \. `# K) g* L9 P7 P
$ ]0 K, F  Q# V

不同阶的导数图像如下:


( j5 z. M1 p' H% T2 b

# G2 i0 V: |5 r" C8 C: j4 g% l' g
定量地分析误差时, 考虑到计算得出的4阶导数向量, 其长度比原始对照向量f4短, 所以两个向量取同样多点进行比较, 就可以得出数值方法的相对误差最大值为

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


9 e8 S5 v/ ]/ F函数实现# @! h2 s8 r  R- p2 X  U
& ?# J: o9 }1 {
  • 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;
    1 y+ U. ~7 H- d3 F) U- Y; l
& ^( I8 {$ w& E0 [% z8 L

  n1 ]0 E- W. r" d1 B, ~! v% ]4 C7 \

该用户从未签到

2#
发表于 2021-1-28 11:10 | 只看该作者
matlab实现数值微分(diff_ctr函数)
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

推荐内容上一条 /1 下一条

EDA365公众号

关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

GMT+8, 2025-10-31 14:35 , Processed in 0.140625 second(s), 26 queries , Gzip On.

深圳市墨知创新科技有限公司

地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

快速回复 返回顶部 返回列表