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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
本帖最后由 piday123 于 2021-1-28 10:25 编辑 2 H5 O: `; v9 V5 e+ |* [

3 Z- u6 Q1 c9 A7 `5 L+ M! E目录
  • 总述
  • 函数说明
  • 应用举例
  • 函数实现+ F. O& ^1 W# e7 T- g" d

* |! @4 |( |) Q3 A; q
& S+ r( {& I3 r5 m# a, e总述/ a. g- D, D- x+ l

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

1 _9 m/ ^) i5 L- N

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


0 C, c9 O# S3 e8 @函数说明
5 f) N* O. K) }+ j4 C0 Z$ U: H! s3 ~: H3 `
  • 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)))
    ' k  a' L  e. A  @  w" I5 i

3 ]& [' U! R/ E* o
8 p' I2 A" ~8 s; c8 i$ p0 D) g
/ _1 b+ x  P3 L4 u& u/ J应用举例

问题: 求函数

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


" W9 \9 H/ @% _) m7 @/ h

代码如下:

% Q( E& y+ V9 h2 K8 P
  • % // 输入函数,并求解析解,并代入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)))" b' \/ V) F" K
+ m& `; T8 }9 ]* P6 j- p

' h0 g6 Z! [& c0 y9 P- y* |/ ~! Y( n. n! z5 B% {

不同阶的导数图像如下:


* e% h; k& E2 A% t( U; Y+ a


1 n% `$ S3 k) n" i: ]8 q2 @: h2 r+ j定量地分析误差时, 考虑到计算得出的4阶导数向量, 其长度比原始对照向量f4短, 所以两个向量取同样多点进行比较, 就可以得出数值方法的相对误差最大值为

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


- d: t' v0 _- _2 Q! ?函数实现1 ~" ?  t7 C" N. F8 w

; F; j4 }" g' C/ f( Q, A" X
  • 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;$ F8 g+ S! @2 M' {- r4 W# H
. i( c2 F- O$ s; e# j

3 p, ~. u+ ]6 p0 z/ e; ~* }: G
  j5 ~5 K2 }: O. l! X' Y4 [

该用户从未签到

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

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-10-6 09:20 , Processed in 0.156250 second(s), 26 queries , Gzip On.

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

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

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