matlab實現數值微分(diff_ctr函式)
阿新 • • 發佈:2020-08-07
目錄
總述
如果已知函式表示式,可以通過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=\frac{sin x}{x^2+4x+3}$ 的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)))
不同階的導數影象如下:
定量地分析誤差時, 考慮到計算得出的4階導數向量, 其長度比原始對照向量f4
短, 所以兩個向量取同樣多點進行比較, 就可以得出數值方法的相對誤差最大值為$3.5 \times 10^{-4}$, 亦即 $0.03 5%$ 。 由此可見, 這裡的數值方法還是很精確的。
函式實現
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;
此函式原始檔可前往下面網址下載: