【文章內(nèi)容簡介】 sin(i*pi*x/L)。 end if a+b, F=subs(F,x,xLa)。 end ％換回變量區(qū)域 例 : syms x。 f=x*(xpi)*(x2*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。 f=abs(x)/x。 % 定義方波信號 xx=[pi:pi/200:pi]。 xx=xx(xx~=0)。 xx=sort([xx,eps,eps])。 % 剔除零點 yy=subs(f,x,xx)。 plot(xx,yy,39。r.39。), hold on % 繪制出理論值并保持坐標系 for n=2:20 [a,b,f1]=fseries(f,x,n), y1=subs(f1,x,xx)。 plot(xx,y1) end a = [ 0, 0, 0] b = [ 4/pi, 0] f1 = 4/pi*sin(x) a = [ 0, 0, 0, 0] b = [ 4/pi, 0, 4/3/pi] f1 = 4/pi*sin(x)+4/3/pi*sin(3*x) …… 級數(shù)求和的計算 ? 是在符號工具箱中提供的 例 :計算 format long。 sum(2.^[0:63]) %數(shù)值計算 ans = +019 sum(sym(2).^[0:200]) % 或 syms k。 symsum(2^k,0,200) ％把 2定義為符號量可使計算更精確 ans = 3213876088517980551083924184682325205044405987565585670602751 syms k。 symsum(2^k,0,200) ans = 3213876088517980551083924184682325205044405987565585670602751 例 :試求解無窮級數(shù)的和 syms n。 s=symsum(1/((3*n2)*(3*n+1)),n,1,inf) %采用符號運算工具箱 s = 1/3 m=1:10000000。 s1=sum(1./((3*m2).*(3*m+1)))。%數(shù)值計算方法 ,雙精度有效位 16,“大數(shù)吃小數(shù)” ,無法精確 format long。 s1 % 以長型方式顯示得出的結(jié)果 s1 = 例 :求解 syms n x s1=symsum(2/((2*n+1)*(2*x+1)^(2*n+1)),n, 0,inf)。 simple(s1) % 對結(jié)果進行化簡， MATLAB 及以前版本因本身 bug 化簡很麻煩 ans = log((((2*x+1)^2)^(1/2)+1)/(((2*x+1)^2)^(1/2)1)) %實際應(yīng)為 log((x+1)/x) 例 :求 syms m n。 limit(symsum(1/m,m,1,n)log(n),n,inf) ans = eulergamma vpa(ans, 70) % 顯示 70 位有效數(shù)字 ans = .5772156649015328606065120900824024310421593359399235988057672348848677 數(shù)值微分 39。0( ) ( ) ( ) l imhf x h f xfxh????差商型求導公式 由導數(shù)定義39。39。39。1 ( ) ( ) ( )2 ( ) ( ) ( )3 ( ) ( ) ( )2f x h f xfxhf x f x hfxhf x h f x hfxh??????? ? ??（）向前差商公式（ ）向后差商公式（）中心差商公式 （中點方法 ） xh x x+h B C A T f(x) 數(shù)值微分算法 ? 向前差商公式： ? 向后差商公式 兩種中心公式： 39。 2 39。39。 3 39。39。39。 439。 2 39。39。 3 39。39。39。 4239。 39。39。39。( ) ( ) ( ) / 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 實現(xiàn) function [dy,dx]=diff_ctr(y, Dt, n) yx1=[y 0 0 0 0 0]。 yx2=[0 y 0 0 0 0]。 yx3=[0 0 y 0 0 0]。 yx4=[0 0 0 y 0 0]。 yx5=[0 0 0 0 y 0]。 yx6=[0 0 0 0 0 y]。 switch n case 1 dy = (diff(yx1)+7*diff(yx2)+7*diff(yx3) … diff(yx4))/(12*Dt)。 L0=3。 case 2 dy=(diff(yx1)+15*diff(yx2) 15*diff(yx3)… +diff(yx4))/(12*Dt^2)。L0=3。 ％數(shù)值計算 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)。 L0=5。 case 4 dy = (diff(yx1)+11*diff(yx2)28*diff(yx3)+28*… diff(yx4)11*diff(yx5)+diff(yx6))/(6*Dt^4)。 L0=5。 end dy=dy(L0+1:endL0)。 dx=([1:length(dy)]+L02(n2))*Dt。 調(diào)用格式： y為 等距實測數(shù)據(jù)， dy為得出的導數(shù)向量， dx為相應(yīng)的自變量向量， dy、 dx的數(shù)據(jù)比 y短 。 [ , ] _ ( , , )yxd d d iff c tr y t n??? 例： 求導數(shù)的解析解 ,再 用數(shù)值微分求取原函數(shù)的 1～4 階導數(shù)，并和解析解比較精度。 h=。 x=0:h:pi。 syms x1。 y=sin(x1)/(x1^2+4*x1+3)。 % 求各階導數(shù)的解析解與對照數(shù)據(jù) yy1=diff(y)。 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=sin(x)./(x.^2+4*x+3)。 % 生成已知數(shù)據(jù)點 [y1,dx1]=diff_ctr(y,h,1)。 subplot(221),plot(x,f1,dx1,y1,39。:39。)。 [y2,dx2]=diff_ctr(y,h,2)。 subplot(222),plot(x,f2,dx2,y2,39。:39。) [y3,dx3]=diff_ctr(y,h,3)。 subplot(223),plot(x,f3,dx3,y3,39。:39。)。 [y4,dx4]=diff_ctr(y,h,4)。 subplot(224),plot(x,f4,dx4,y4,39。:39。) 求最大相對誤差： norm((y4… f4(4:60))./f4(4:60)) ans = 用插值、擬合多項式的求導數(shù) ? 基本思想：當已知函數(shù)在一些離散點上的函數(shù)值時，該函數(shù)可用插值或擬合多項式來近似，然后對多項式進行微分求得導數(shù)。 ? 選取 x=0附近的少量點 ? 進行多項式擬合或插值 ? g(x)在 x=0處的 k階導數(shù)為 ( , ) , 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????? 通過坐標變換用上述方法計算任意 x點處的導數(shù)值 ? 令 ? 將 g(x)寫成 z的表達式 ? 導數(shù)為 ? 可直接用 擬合節(jié)點 得到系數(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: 計算 x=a=。 xd=[ 0 ]。 yd=[ ]。 a=。L=length(xd)。 d=polyfit(xda,yd,L1)。fact=[1]。 for k