关于MATLAB数值微分

yxzjs

      对于离散数值，在MATLAB中没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其调用格式为：

DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,…,n-1。

DX=diff(X,n)：计算X的n阶向前差分。例如，diff(X,2)=diff(diff(X))。

DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。


      还有一个gradient函数可以直接求矩阵的二维差分：
[dy,dx]=gradient(A)


      实际上无论是向前差分还是向后差分的数值算法，在求解高阶微分时，其精度是较低的。不知在MATLAB中，有较高精度的中心差分算法有没有现成的函数可用？

bainhome

应该可以尝试使用样条求导命令fnder计算，具体请查帮助。但是它的算法不是中心差分，样条类命令的算法好像要复杂一点儿。
中心差分算法在薛定宇的《高等应用数学问题的MATLAB求解》书中有代码(四阶代数精度)：

1.    
   
2. function [dy,dx]=diff_ctr(y, Dt, n)
   
3. yx1=[y 0 0 0 0 0]; yx2=[0 y 0 0 0 0]; yx3=[0 0 y 0 0 0];
   
4. yx4=[0 0 0 y 0 0]; yx5=[0 0 0 0 y 0]; yx6=[0 0 0 0 0 y];
   
5. switch n
   
6.       case 1
   
7.          dy = (-diff(yx1)+7*diff(yx2)+7*diff(yx3)-diff(yx4))/(12*Dt); L0=3;
   
8.       case 2
   
9.          dy=(-diff(yx1)+15*diff(yx2)-15*diff(yx3)+diff(yx4))/(12*Dt^2);L0=3;
   
10.       case 3
    
11.          dy=(-diff(yx1)+7*diff(yx2)-6*diff(yx3)-6*diff(yx4)+...
    
12.                7*diff(yx5)-diff(yx6))/(8*Dt^3); L0=5;
    
13.       case 4
    
14.          dy = (-diff(yx1)+11*diff(yx2)-28*diff(yx3)+28*diff(yx4)-...
    
15.             11*diff(yx5)+diff(yx6))/(6*Dt^4);L0=5;
    
16. end
    
17. dy=dy(L0+1:end-L0); dx=([1:length(dy)]+L0-2-(n>2))*Dt;
    
18. dy=dy(L0+1:end-L0); dx=[0:length(dy)-1]*Dt;

复制代码

测试代码：

1. h=0.05; x=0:h:pi; syms x1; y=sin(x1)/(x1^2+4*x1+3);
   
2. yy1=diff(y); f1=subs(yy1,x1,x);   % 求各阶导数的解析解与对照数据
   
3. yy2=diff(yy1); f2=subs(yy2,x1,x); yy3=diff(yy2); f3=subs(yy3,x1,x);
   
4. yy4=diff(yy3); f4=subs(yy4,x1,x);
   
5. y=sin(x)./(x.^2+4*x+3);   % 生成已知数据点
   
6. [y1,dx1]=diff_ctr(y,h,1); subplot(221),plot(x,f1,dx1,y1,':');
   
7. [y2,dx2]=diff_ctr(y,h,2); subplot(222),plot(x,f2,dx2,y2,':')
   
8. [y3,dx3]=diff_ctr(y,h,3); subplot(223),plot(x,f3,dx3,y3,':');
   
9. [y4,dx4]=diff_ctr(y,h,4); subplot(224),plot(x,f4,dx4,y4,':')
   
10. norm((y4-f4(4:60))./f4(4:60))

复制代码

至于gradient命令，记得好像以前见过一个解释的例子，应该也是在薛定宇的书中，说明其结果并不是真正的微分结果，还要具体除以网格点步距。

[ 本帖最后由 bainhome 于 2008-1-27 02:01 编辑 ]

yxzjs

多谢bainhome提供的中心差分算法。同时也提醒了我：使用diff命令后也得除以网格点步距，难怪我的结果对不上。

sigma665

你的diff_ctr比书上多了最后一行，这一行有什么作用呢