ch2 - 符号计算2010a(4)

6 21 32.82.62.42.221.81.61.41.2101234567

图2.2-2

12

2.3

2.3.1

符号微积分

极限和导数的符号计算

【例2.3-1】

syms t x k

s=sin(k*t)/(k*t); f=(1-1/x)^(k*x);

Lsk=limit(s,0) Ls1=subs(Lsk,k,1) Lf=limit(f,x,inf)

Lf1=vpa(subs(Lf,k,sym('-1')),48) Lsk = 1

Ls1 = 1 Lf =

1/exp(k) Lf1 =

2.7182818284590452353602874713526624977572470937

【例2.3-2】

syms a t x;

f=[a,t^3;t*cos(x), log(x)]; df=diff(f) dfdt2=diff(f,t,2) dfdxdt=diff(diff(f,x),t) df =

[ 0, 0] [ -t*sin(x), 1/x] dfdt2 = [ 0, 6*t] [ 0, 0] dfdxdt =

[ 0, 0] [ -sin(x), 0]

【例2.3-3】

syms x1 x2;

f=[x1*exp(x2);x2;cos(x1)*sin(x2)]; v=[x1;x2]; Jf=jacobian(f,v) Jf =

[ exp(x2), x1*exp(x2)] [ 0, 1] [ -sin(x1)*sin(x2), cos(x1)*cos(x2)]

【例2.3-4】 （1）

clear syms x

syms d positive

f_p=sin(x); df_p=limit((subs(f_p,x,x+d)-f_p)/d,d,0) % <5> df_p0=limit((subs(f_p,x,d)-subs(f_p,x,0))/d,d,0)% <6> df_p = cos(x) df_p0 =

13

1

（2）

f_n=sin(-x); df_n=limit((f_n-subs(f_n,x,x-d))/d,d,0) % <8> df_n0=limit((subs(f_n,x,0)-subs(f_n,x,-d))/d,d,0) % <9> df_n = -cos(x) df_n0 = -1

（3）

f=sin(abs(x));

dfdx=diff(f,x) dfdx0=subs(dfdx,x,0) dfdx =

cos(abs(x))*sign(x) dfdx0 = 0

% <12>

% <11>

（4）

clf

xn=-3/2*pi:pi/50:0;xp=0:pi/50:3/2*pi;xnp=[xn,xp(2:end)]; hold on

plot(xnp,subs(f,x,xnp),'k','LineWidth',3) plot(xn,subs(df_n,x,xn),'--r','LineWidth',3) plot(xp,subs(df_p,x,xp),':r','LineWidth',3)

legend(char(f),char(df_n),char(df_p),'Location','NorthEast') grid on

xlabel('x') hold off 10.80.60.40.20-0.2-0.4-0.6-0.8-1 -5-4-3-2-10x12345sin(abs(x))-cos(x)cos(x) % <14> % <17>

图 2.3-1

【例2.3-5】 （1）

14

clear syms x

g=sym('cos(x+sin(y(x)))=sin(y(x))') % <3> dgdx=diff(g,x) % <4> g =

cos(x + sin(y(x))) = sin(y(x)) dgdx = -sin(x + sin(y(x)))*(cos(y(x))*diff(y(x), x) + 1) = cos(y(x))*diff(y(x), x)

（2）

dgdx1=subs(dgdx,'diff(y(x),x)','dydx') % <5> dgdx1 =

-sin(x + sin(y(x)))*(dydx*cos(y(x)) + 1) = dydx*cos(y(x))

（3）

dydx=solve(dgdx1,'dydx') dydx =

-sin(x + sin(y(x)))/(cos(y(x)) + cos(y(x))*sin(x + sin(y(x))))

【例2.3-6】 （1）

syms x

r=taylor(x*exp(x),9,x,0) % <2> pretty(r) % r =

x^8/5040 + x^7/720 + x^6/120 + x^5/24 + x^4/6 + x^3/2 + x^2 + x

8 7 6 5 4 3

x x x x x x 2

---- + --- + --- + -- + -- + -- + x + x 5040 720 120 24 6 2

（2）

R = evalin(symengine,'series(x*exp(x),x=0,8)') % <4> pretty(R) % R =

x + x^2 + x^3/2 + x^4/6 + x^5/24 + x^6/120 + x^7/720 + x^8/5040 + O(x^9)

3 4 5 6 7 8

2 x x x x x x 9 x + x + -- + -- + -- + --- + --- + ---- + O(x ) 2 6 24 120 720 5040

【例2.3-7】

TL1=evalin(symengine,'mtaylor(sin(x^2+y),[x,y],8)') TL1 =

(x^6*y^2)/12 - x^6/6 + (x^4*y^3)/12 - (x^4*y)/2 - (x^2*y^6)/720 + (x^2*y^4)/24 - (x^2*y^2)/2 + x^2 - y^7/5040 + y^5/120 - y^3/6 + y

class(TL1) ans = sym

2.3.2

【例2.3-8】。

序列/级数的符号求和

15