【正文】 ( , ) , 1 , 2 , , 1iix y i n??11 2 1()nnnng x c x c x c x c??? ? ? ? ?() 1( 0 ) ! 0 , 1 , 2 , ,k nkg c k k n????? 通過坐標(biāo)變換用上述方法計(jì)算任意 x點(diǎn)處的導(dǎo)數(shù)值 ? 令 ? 將 g(x)寫成 z的表達(dá)式 ? 導(dǎo)數(shù)為 ? 可直接用 擬合節(jié)點(diǎn) 得到系數(shù) d=polyfit(xa,y,length(xd)1) z x a??11 2 1( ) ( )nnnng x g z d z d z d z d??? ? ? ? ? ?( ) ( )1( ) ( 0 ) ! 0 , 1 , ,kknkg a g d k k n??? ? ?()gz ( , )iix a y? id? 例：數(shù)據(jù)集合如下： xd: 0 yd: 計(jì)算 x=a=。:39。 [y4,dx4]=diff_ctr(y,h,4)。:39。) [y3,dx3]=diff_ctr(y,h,3)。 subplot(222),plot(x,f2,dx2,y2,39。)。 subplot(221),plot(x,f1,dx1,y1,39。 y=sin(x)./(x.^2+4*x+3)。 yy4=diff(yy3)。 yy3=diff(yy2)。 yy2=diff(yy1)。 % 求各階導(dǎo)數(shù)的解析解與對(duì)照數(shù)據(jù) yy1=diff(y)。 syms x1。 h=。 調(diào)用格式： y為 等距實(shí)測(cè)數(shù)據(jù)， dy為得出的導(dǎo)數(shù)向量， dx為相應(yīng)的自變量向量， dy、 dx的數(shù)據(jù)比 y短 。 end dy=dy(L0+1:endL0)。 case 4 dy = (diff(yx1)+11*diff(yx2)28*diff(yx3)+28*… diff(yx4)11*diff(yx5)+diff(yx6))/(6*Dt^4)。 ％數(shù)值計(jì)算 diff(X)表示數(shù)組 X相鄰兩數(shù)的差 case 3 dy=(diff(yx1)+7*diff(yx2)6*diff(yx3)6*diff(yx4)+... 7*diff(yx5)diff(yx6))/(8*Dt^3)。 case 2 dy=(diff(yx1)+15*diff(yx2) 15*diff(yx3)… +diff(yx4))/(12*Dt^2)。 switch n case 1 dy = (diff(yx1)+7*diff(yx2)+7*diff(yx3) … diff(yx4))/(12*Dt)。 yx5=[0 0 0 0 y 0]。 yx3=[0 0 y 0 0 0]。( ) ( ) ( ) / 2 ! ( ) / 3! ( )()2( ) ( ) ( ) / 2 ! ( ) / 3! ( )2( ) ( )3!f x t f x t f x t f o tfxtf x t f x t f x t f o tttf x f???? ? ? ? ? ? ? ???? ? ? ? ? ? ? ?????? 中心差分方法及其 MATLAB 實(shí)現(xiàn) function [dy,dx]=diff_ctr(y, Dt, n) yx1=[y 0 0 0 0 0]。39。 4239。39。39。 439。39。39。1 ( ) ( ) ( )2 ( ) ( ) ( )3 ( ) ( ) ( )2f x h f xfxhf x f x hfxhf x h f x hfxh??????? ? ??（）向前差商公式（ ）向后差商公式（）中心差商公式 （中點(diǎn)方法 ） xh x x+h B C A T f(x) 數(shù)值微分算法 ? 向前差商公式： ? 向后差商公式 兩種中心公式： 39。39。 limit(symsum(1/m,m,1,n)log(n),n,inf) ans = eulergamma vpa(ans, 70) % 顯示 70 位有效數(shù)字 ans = .5772156649015328606065120900824024310421593359399235988057672348848677 數(shù)值微分 39。 s1 % 以長(zhǎng)型方式顯示得出的結(jié)果 s1 = 例 :求解 syms n x s1=symsum(2/((2*n+1)*(2*x+1)^(2*n+1)),n, 0,inf)。 s1=sum(1./((3*m2).*(3*m+1)))。 symsum(2^k,0,200) ans = 3213876088517980551083924184682325205044405987565585670602751 例 :試求解無窮級(jí)數(shù)的和 syms n。 sum(2.^[0:63]) %數(shù)值計(jì)算 ans = +019 sum(sym(2).^[0:200]) % 或 syms k。), hold on % 繪制出理論值并保持坐標(biāo)系 for n=2:20 [a,b,f1]=fseries(f,x,n), y1=subs(f1,x,xx)。 plot(xx,yy,39。 xx=sort([xx,eps,eps])。 % 定義方波信號(hào) xx=[pi:pi/200:pi]。 [A,B,F]=fseries(f,x,6,0,2*pi) A = [ 0, 0, 0, 0, 0, 0, 0] B = [ 12, 3/2, 4/9, 3/16, 12/125, 1/18] F = 12*sin(x)+3/2*sin(2*x)+4/9*sin(3*x)+3/16*sin(4*x)+12/125*sin(5*x)+1/18*sin(6*x) 例 : syms x。 end ％換回變量區(qū)域 例 : syms x。 F=F+an*cos(i*pi*x/L)+bn*sin(i*pi*x/L)。 A=[A, an]。 %計(jì)算 a0 for i=1:n an=int(f*cos(i*pi*x/L),x,L,L)/L。 B=[]。 if a+b, f=subs(f,x,x+L+a)。 b=pi。[x=a]39。mtaylor39。,3)。,f,39。 F=maple(39。[x,y]39。mtaylor39。,8) F = mtaylor((x^22*x)*exp(x^2y^2x*y),[x, y],8) maple(‘readlib(mtaylor)’)。,f,39。 F=maple(39。 line(x0,y1) end p = x1/6*x^3+1/120*x^51/5040*x^7 p = x1/6*x^3+1/120*x^51/5040*x^7+1/362880*x^9 p = x1/6*x^3+1/120*x^51/5040*x^7+1/362880*x^91/39916800*x^11 p = x1/6*x^3+1/120*x^51/5040*x^7+1/362880*x^91/39916800*x^11+1/6227020800*x^13 p = x1/6*x^3+1/120*x^51/5040*x^7+1/362880*x^91/39916800*x^11+1/6227020800*x^131/1307674368000*x^15 多變量函數(shù)的 Taylor 冪級(jí)數(shù)展開 ? 多變量函數(shù) 在 的 Taylor冪級(jí)數(shù)的展開 12( , , , )nf x x x 12( , , , )na a a? 例： ？？？ syms x y。), axis([2*pi,2*pi,])。 plot(x0,y0,39。 syms x。 x0=2*pi::2*pi。 pretty(y1) 2 23 3 34 4 4087 5 3067 6 515273 7 386459 8 1/3 x 4/9 x + x x + x x + x x 54 81 9720 7290 1224720 918540 taylor(f,x,9,2) ans = 1/15*sin(2)+(1/15*cos(2)8/225*sin(2))*(x2)+ (127/6750*sin(2)8/225*cos(2))*(x2)^2 +(23/6750*cos(2)+628/50625*sin(2))*(x2)^3 +(15697/6075000*sin(2)+28/50625*cos(2))*(x2)^4 +(203/6075000*cos(2)+6277/11390625*sin(2))*(x2)^5 +(585671/2733750000*sin(2)623/11390625*cos(2))*(x2)^6 +(262453/19136250000*cos(2)+397361/5125781250*sin(2))*(x2)^7 +(875225059/34445250000000*sin(2)131623/35880468750*cos(2))*(x2)^8 syms a。 f=sin(x)/(x^2+4*x+3)。這是由于