ch2 - 符号计算2010a

第 2 章 符号计算

2.1

2.1.1 1 2

符号对象和符号表达式

符号对象的创建和衍生 生成符号对象的基本规则 符号数字和符号常数

【例2.1-1】

a=pi+sqrt(5)

sa=sym('pi+sqrt(5)') Ca=class(a) Csa=class(sa) vpa(sa-a) a =

5.3777 sa =

pi + 5^(1/2) Ca = double Csa = sym ans =

0.000000000000000013822375841085200048593542564188

3 4

基本符号变量 自由符号变量

【例2.1-2】。 （1）

syms u v w z a5 f=sym('3');

Eq=sin(f)*u*z^2+v*z+f*w-a5;

（2）

symvar(Eq) ans =

[ a5, u, v, w, z]

symvar(Eq,100) ans =

[ w, z, v, u, a5]

symvar(Eq,1) ans = w

（3）

result_1=solve(Eq) result_1 =

a5/3 - (v*z)/3 - (u*sin(3)*z^2)/3

1

（4）

result_2=solve(Eq,z) result_2 =

-(v - (v^2 + 4*a5*u*sin(3) - 12*u*w*sin(3))^(1/2))/(2*u*sin(3)) -(v + (v^2 + 4*a5*u*sin(3) - 12*u*w*sin(3))^(1/2))/(2*u*sin(3))

【例2.1-3】 （1）

syms a b x X Y k=sym('3'); z=sym('c*sqrt(d)+y*sin(t)'); EXPR=a*z*X+(b*x^2+k)*Y;

（2）

symvar(EXPR) ans =

[ X, Y, a, b, c, d, t, x, y]

（3）

symvar(EXPR,10) ans =

[ x, y, t, d, c, b, a, X, Y]

（4）

symvar(EXPR,1) ans = x

（5）

symvar(EXPR,3) ans =

[ x, y, t]

（6）

E3=sym('a*sqrt(theta)')

??? Error using ==> sym.sym>sym.sym/scalarsym at 382 Error: argument must be of 'Type::Arithmetical' [sqrt]

Error in ==> sym.sym>sym.sym/char2sym at 337 Scell = scalarsym(x);

Error in ==> sym.sym>sym.sym/symchar at 182 Scell = char2sym(x);

Error in ==> sym.sym>sym.sym at 120

S = cell2sym(S,symchar(x,a,nargin));

E4=sym('a*sqrt(theta123)') E4 =

a*theta123^(1/2)

E5=sym('a*sqrt(theta*t)') E5 =

a*(t*theta)^(1/2)

【例2.1-4】

syms a b t u v x y

2

A=[a+b*x,sin(t)+u;x*exp(-t),log(y)+v] symvar(A,1) A =

[ a + b*x, u + sin(t)] [ x/exp(t), v + log(y)] ans = x

2.1.2 2.1.3 2.1.4

【例2.1-5】。 （1）

符号计算中的算符 符号计算中的函数指令 符号对象的识别

clear

a=1;b=2;c=3;d=4; Mn=[a,b;c,d] Mc='[a,b;c,d]' Ms=sym(Mc) Mn =

1 2 3 4 Mc =

[a,b;c,d] Ms = [ a, b] [ c, d]

（2）

SizeMn=size(Mn) SizeMc=size(Mc) SizeMs=size(Ms) SizeMn =

2 2 SizeMc =

1 9 SizeMs =

2 2

（3）

CMn=class(Mn) CMc=class(Mc) CMs=class(Ms) CMn = double CMc = char CMs = sym

（4）

isa(Mn,'double') isa(Mc,'char') isa(Ms,'sym') ans = 1 ans =

3

1 ans = 1

（5）

whos Mn Mc Ms

Name Size Bytes Class Attributes

Mc 1x9 18 char Mn 2x2 32 double Ms 2x2 60 sym

2.1.5 1 2 3

符号运算机理和变量假设 符号运算的工作机理 对符号变量的限定性假设 清除变量和撤销假设

【例2.1-6】 （1）

syms x clear

f=x^3+4.75*x+2.5; rf=solve(f,x) rf =

-1/2 1/4 - (79^(1/2)*i)/4 (79^(1/2)*i)/4 + 1/4

evalin(symengine,'getprop(x)') ans = C_

（2）

syms x real

rfr=solve(f,x) rfr = -1/2

evalin(symengine,'getprop(x)') ans = R_

（3）

clear x syms x

g=x^2+x+5; rg=solve(g,x)

Warning: Explicit solution could not be found. > In solve at 98 rg =

[ empty sym ]

（4）

syms x clear rg=solve(g,x) rg =

4

- (19^(1/2)*i)/2 - 1/2 (19^(1/2)*i)/2 - 1/2

【例2.1-7】 （1）

clear all reset(symengine) Da=1.2;Dw=1/3; syms sa sw sx sy sz syms A B positive syms C real

（2）

whos Name Size Bytes Class Attributes

A 1x1 58 sym B 1x1 58 sym C 1x1 58 sym Da 1x1 8 double Dw 1x1 8 double sa 1x1 60 sym sw 1x1 60 sym sx 1x1 60 sym sy 1x1 60 sym sz 1x1 60 sym

（3）

syms

'A' 'B' 'C' 'sa' 'sw' 'sx' 'sy' 'sz'

（4）

evalin(symengine,'anames(Properties)')

ans =

{A, B, C}

（5）

clear A syms

'B' 'C' 'ans' 'sa' 'sw' 'sx' 'sy' 'sz'

evalin(symengine,'anames(Properties)') ans =

{A, B, C}

（6）

syms B clear syms

'B' 'C' 'ans' 'sa' 'sw' 'sx' 'sy' 'sz'

evalin(symengine,'anames(Properties)') ans = {A, C}

5

2.1.6

【例2.1-8】 （1）

符号帮助体系

图2.1-1

（2）

6

图2.1-2

? （3）

7

图2.1-3

2.2

符号数字及表达式的操作

双精度数字与符号数字之间的转换 双精度数字向符号数字的转换 符号数字向双精度数字转换

符号数字的任意精度表达形式

2.2.1 1 2 2.2.2

【例2.2-1】 （1）

reset(symengine)

sa=sym('1/3+sqrt(2)') sa =

2^(1/2) + 1/3

（2）

digits

8

Digits = 32

format long a=1/3+sqrt(2) sa_Plus_a=vpa(sa+a,20) sa_Minus_a=vpa(sa-a,20) a =

1.747546895706428 sa_Plus_a =

3.4950937914128567869 sa_Minus_a =

-0.000000000000000022658064826339973669

（3）

sa32=vpa(sa) % digits(48) % <9> sa5=vpa(sa,5) % <10>

sa48=vpa(sa) % <11> sa32 =

1.747546895706428382135022057543 sa5 = 1.7475 sa48 =

1.74754689570642838213502205754303141190300520871

（4）

a

a2=a+1/10^12 a_16=vpa(a,16);a2_16=vpa(a2,16); a_12=vpa(a,12);a2_12=vpa(a2,12); a =

1.747546895706428 a2 =

1.747546895707429

a_16==a2_16 ans = 0

a_12==a2_12 ans = 0

2.2.3

符号表达式的基本操作

【例2.2-2】

syms x

f=(1/x^3+6/x^2+12/x+8)^(1/3) g1=simple(f) f =

(12/x + 6/x^2 + 1/x^3 + 8)^(1/3) g1 =

((2*x + 1)^3/x^3)^(1/3)

9

<8>

2.2.4 1

表达式中的置换操作 公因子法简化表达

【例2.2-3】。 （１）

clear A=sym('[a b;c d]') [V,D]=eig(A) A =

[ a, b] [ c, d] V =

[(a/2+d/2-a^2-2*a*d+d^2+4*b*c)^(1/2)/2)/c-d/c, (a/2+d/2+(a^2-2*a*d+ d^2+4*b*c)^(1/2)/2)/c-d/c] [ 1, 1] D = [a/2 + d/2 - (a^2 - 2*a*d + d^2 + 4*b*c)^(1/2)/2, 0] [ 0, a/2 + d/2 + (a^2 - 2*a*d + d^2 + 4*b*c)^(1/2)/2]

（2）

subexpr([V;D]) % <4> who % sigma =

(a^2 - 2*a*d + d^2 + 4*b*c)^(1/2) ans =

[ (a/2 + d/2 - sigma/2)/c - d/c, (a/2 + d/2 + sigma/2)/c - d/c] [ 1, 1] [ a/2 + d/2 - sigma/2, 0] [ 0, a/2 + d/2 + sigma/2]

Your variables are:

A D V ans sigma

（3）

Dw=subexpr(D,'w') w =

(a^2 - 2*a*d + d^2 + 4*b*c)^(1/2) Dw =

[ a/2 + d/2 - w/2, 0] [ 0, a/2 + d/2 + w/2]

（4）

[RVD,w]=subexpr([V;D],'w') % <7> RVD =

[ (a/2 + d/2 - w/2)/c - d/c, (a/2 + d/2 + w/2)/c - d/c] [ 1, 1] [ a/2 + d/2 - w/2, 0] [ 0, a/2 + d/2 + w/2] w =

(a^2 - 2*a*d + d^2 + 4*b*c)^(1/2)

2 通用置换指令

【例2.2-4】 （1）

clear

10

syms a b x; f=a*sin(x)+b f =

b + a*sin(x)

（2）

f1=subs(f,sin(x),'log(y)') class(f1)

f1 =

b + a*log(y) ans = sym

（3）

f2=subs(f,a,3.11) class(f2)

f2 =

b + (311*sin(x))/100 ans = sym

（4）

f3=subs(f,{a,b,x},{2,5,sym('pi/3')})

class(f3) f3 =

3^(1/2) + 5 ans = sym

（5）

format format compact f4=subs(f,{a,b,x},{2,5,pi/3})

class(f4) f4 =

6.7321 ans = double

（6）

f5=subs(f,x,0:pi/2:pi) class(f5) f5 =

[ b, a + b, b] ans = sym

（7）

t=0:pi/10:2*pi;

f6=subs(f,{a,b,x},{2,3,t}) plot(t,f6) f6 =

Columns 1 through 8

3.0000 3.6180 4.1756 4.6180 4.9021 5.0000 4.9021 4.6180

Columns 9 through 16

4.1756 3.6180 3.0000 2.3820 1.8244 1.3820 1.0979 1.0000

Columns 17 through 21

1.0979 1.3820 1.8244 2.3820 3.0000

11

54.543.532.521.5101234567图2.2-1

（8）

k=(0.5:0.1:1)';

f6=subs(subs(f,{a,b},{k,2}),x,t); size(f6) plot(t,f6) ans =

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

（1）

syms n k

f1=1/(k*(k+1));

s1=symsum(f1,k,1,n) s1 =

1 - 1/(n + 1)

（2）

f2=x^(2*k-1)/(2*k-1); s2=symsum(f2,k,1,inf) s2 =

piecewise([abs(x) < 1, atanh(x)])

（3）

f3=[1/(2*k-1)^2,(-1)^k/k]; s3=symsum(f3,k,1,inf) s3 =

[ pi^2/8, -log(2)]

2.3.3

符号积分

【例2.3-9】

clear

syms a b x f1=x*log(x) s1=int(f1,x) s1=simple(s1) % f1 =

x*log(x) s1 =

(x^2*(log(x) - 1/2))/2 s1 =

x^2*(log(x)/2 - 1/4)

【例2.3-10】

f2=[a*x,b*x^2;1/x,sin(x)] disp(' ')

disp('The integral of f is') pretty(int(f2)) f2 =

[ a*x, b*x^2] [ 1/x, sin(x)]

The integral of f is

+- -+ | 2 3 | | a x b x | | ----, ---- | | 2 3 | | | | log(x), -cos(x) | +- -+

【例2.3-11】

syms x y z

16

<5>

F2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2) VF2=vpa(F2)

Warning: Explicit integral could not be found. F2 =

(14912*2^(1/4))/4641 - (6072064*2^(1/2))/348075 + (64*2^(3/4))/225 + 1610027357/6563700 VF2 =

224.921535733311431597907100328046757677071376012

【例2.3-12】。 （1）

syms a r theta phi

r=a*theta; x=r*cos(theta); y=r*sin(theta);

dLdth=sqrt(diff(x,theta)^2+diff(y,theta)^2); warning off % <5>

L=simple(int(dLdth,theta,0,phi)) L =

(phi*(a^2*phi^2 + a^2)^(1/2) + log(phi +

1)^(1/2))*(a^2)^(1/2))/2

（2）

L_2pi=subs(L,[a,phi],sym('[1,2*pi]')) L_2pi_vpa=vpa(L_2pi) L_2pi =

log(2*pi + (4*pi^2 + 1)^(1/2))/2 + pi*(4*pi^2 + 1)^(1/2) L_2pi_vpa =

21.2562941482090988007025122725661088234709310476

（3）

L1=subs(L,a,sym('1'));

ezplot(L1*cos(phi),L1*sin(phi),[0,2*pi]) grid on hold on x1=subs(x,a,sym('1')); y1=subs(y,a,sym('1')); h1=ezplot(x1,y1,[0,2*pi]); % <14> set(h1,'Color','r','LineWidth',5) %

title(' ') % <16> legend('螺线长度-幅角曲线','阿基米德螺线') hold off %

17

% <6>

+

(phi^2

420-2-4y螺线长度-幅角曲线阿基米德螺线-6-8-10-12-14-16 -505x101520图 2.3-2 阿基米德螺线（粗红）和螺线长度函数（细蓝）

2.4

2.4.1 2.4.2

微分方程的符号解法

符号解法和数值解法的互补作用 求微分方程符号解的一般指令 微分方程符号解示例

2.4.3

【例2.4-1】

clear all % <1> S=dsolve('Dx=y,Dy=-x') % disp(' ')

disp(['微分方程的解',blanks(2),'x',blanks(22),'y']) disp([S.x,S.y]) S =

x: [1x1 sym] y: [1x1 sym]

微分方程的解 x y

[ C2*cos(t) + C1*sin(t), C1*cos(t) - C2*sin(t)]

【例2.4-2】 （1）

clear all % <1>

y=dsolve('(Dy)^2-x*Dy+y=0','x') % <2> y =

18

x^2/4 C3*x - C3^2

（2）

clf,hold on %

hy1=ezplot(y(1),[-6,6,-4,8],1); % <4> set(hy1,'Color','r','LineWidth',5)

for k=-2:0.5:2 % <6> y2=subs(y(2),'C3',k); %<7> ezplot(y2,[-6,6,-4,8],1) end % <9> hold off % box on

legend('奇解','通解','Location','Best') ylabel('y')

title(['\\fontsize{14}微分方程',' (y '')^2 – xy '' + y = 0 ','的解'])

微分方程 (y ')2 – xy ' + y = 0 的解8奇解通解6 42y0-2-4 -6-4-20x246图 2.4-1 通解和奇解曲线

【例2.4-3】 （1）

y=dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=0','x') y =

(31*x^4)/468 - x^3/3 + 125/468

（2）

xn=-1:6;

yn=subs(y,'x',xn) ezplot(y,[-1,6]) hold on

plot([1,5],[0,0],'.r','MarkerSize',20) text(1,1,'y(1)=0') text(4,1,'y(5)=0')

title(['x*D2y - 3*Dy = x^2',', y(1)=0,y(5)=0']) hold off yn =

0.6667 0.2671 0 -1.3397 -3.3675 -4.1090 0.0000

19

14.1132 x*D2y - 3*Dy = x2, y(1)=0,y(5)=0543210-1-2-3-4-1012x3456y(1)=0y(5)=0图 2.4-2 两点边值问题的解曲线

2.5

2.5.1

符号变换和符号卷积

Fourier变换及其反变换

【例 2.5-1】 （1）

syms t w

ut=heaviside(t); UT=fourier(ut) UT =

pi*dirac(-w) - i/w

（2）

Ut=ifourier(UT,w,t)

SUt=simple(Ut) Ut =

(pi + pi*(2*heaviside(t) - 1))/(2*pi) SUt =

heaviside(t)

%

<5>

（3）

t=-2:0.01:2;

ut=heaviside(t);

kk=find(t==0); % <8> plot(t(kk),ut(kk),'.r','MarkerSize',30) hold on

ut(kk)=NaN; % <10> plot(t,ut,'-r','LineWidth',3)

plot([t(kk),t(kk)],[ut(kk-1),ut(kk+1)],'or','MarkerSize',10) %

20

hold off grid on

axis([-2,2,-0.2,1.2])

xlabel('\\fontsize{14}t'),ylabel('\\fontsize{14}ut') title('\\fontsize{14}Heaviside(t)') Heaviside(t)10.80.6ut0.40.20-0.2-2-1.5-1-0.500.511.52t图 2.5-1 Heaviside(t) 定义的单位阶跃函数

【例2.5-2】 （1）

syms A t w tao

yt=A*(heaviside(t+tao/2)-heaviside(t-tao/2)); % Yw=fourier(yt,t,w) % Yw_fy=simplify(Yw) % Yw_fy_e=simple(Yw_fy) % Yw =

A*((1/exp((tao*w*i)/2))*(- pi*dirac(-w) + i/w) - exp((tao*w*i)/2)*(- pi*dirac(-w) + i/w)) Yw_fy =

(2*A*sin((tao*w)/2))/w Yw_fy_e =

(2*A*sin((tao*w)/2))/w

（2）

Yt=ifourier(Yw_fy_e,w,t) %

Yt_e=simple(Yt) % <7> Yt =

(A*transform::fourier(sin((tao*w)/2)/w, w, t))/pi Yt_e =

-A*(heaviside(t - tao/2) - heaviside(t + 1/2*tao))

（3）

t3=3;

tn=-3:0.1:3;

yt13=subs(yt,{A,tao},{1,t3}) yt13n=subs(yt13,'t',tn);

21

kk=find(tn==-t3/2|tn==t3/2); % <13>

plot(tn(kk),yt13n(kk),'.r','MarkerSize',30) yt13n(kk)=NaN; % <15> hold on

plot(tn,yt13n,'-r','LineWidth',3) hold off grid on

axis([-3,3,-0.5,1.5]) yt13 =

heaviside(t + 3/2) - heaviside(t - 3/2) 1.5%

10.50-0.5-3-2-10123图 2.5-2 由 Heaviside(t) 构造的矩形波

（4）

Yw13=subs(Yw_fy_e,{A,tao},{1,t3}); subplot(2,1,1),ezplot(Yw13),grid on

subplot(2,1,2),ezplot(abs(Yw13)),grid on

22

(2 sin((3 w)/2))/w3210-1-6-402w(2 abs(sin((3 w)/2)))/abs(w)-2463210-6-4-20w246图 2.5-3 矩形脉冲的频率曲线和幅度频谱

【例2.5-3】 （1）

clear

syms t x w

ft=exp(-(t-x))*heaviside(t-x); gt=exp(-(t-x));

% %

<3> <4>

（2）

F1=simple(fourier(ft,t,w)) % G1=simple(fourier(gt,t,w)) % F1 =

(1/exp(w*x*i))/(1 + w*i) G1 =

transform::fourier(exp(x - t), t, -w)

（3）

F2=simple(fourier(ft,t)) % F3=simple(fourier(ft)) % F2 =

-exp(t^2*i)/(- 1 + t*i) F3 =

-(1/exp(t*w*i))/(- 1 + w*i)

2.5.2

【例2.5-4】 （1）

Laplace变换及其反变换

syms t s a b

f1=exp(-a*t)*sin(b*t) % <2>

23

F1=laplace(f1,t,s) f1 =

sin(b*t)/exp(a*t) F1 =

b/((a + s)^2 + b^2)

（2）

sym a clear %

<4>

f2=heaviside(t-a)

F2=laplace(f2,t,s) % ans = a f2 =

heaviside(t - a) F2 =

laplace(heaviside(t - a), t, s)

syms a positive % F3=laplace(f2) %

F3 =

1/(s*exp(a*s))

（4）

f4=dirac(t-b); % F4=laplace(f4,t,s) % F4 =

piecewise([b < 0, 0], [0 <= b, 1/exp(b*s)])

f5=dirac(t-a); % <11> F5=laplace(f5,t,s) % ft_F5=ilaplace(F5,s,t) % F5 =

1/exp(a*s) ft_F5 =

dirac(a - t)

（5）

n=sym('n','clear'); % <14> F6=laplace(t^n,t,s) % F6 =

piecewise([-1 < Re(n), gamma(n + 1)/s^(n + 1)])

n=sym('n','positive') % <16>

F6=laplace(t^n,t,s) % <17> n = n F6 =

gamma(n + 1)/s^(n + 1)

2.5.3

Z变换及其反变换

【例2.5-5】 （1）

clear

syms n z clear % <2>

gn=6*(1-(1/2)^n) % G=simple(ztrans(gn,n,z)); % pretty(G) %

gn =

6 - 6*(1/2)^n

24

<7>

6 z

-------------- 2

2 z - 3 z + 1

（2）

syms n w T z clear fwn=sin(w*n*T); FW=ztrans(fwn,n,z); pretty(FW),disp(' ') inv_FW=iztrans(FW,z,n)

z sin(T w)

--------------------- 2

z - 2 cos(T w) z + 1

inv_FW =

sin(T*n*w)

% % % %

<6>

（3）

syms n z clear f1=1;

F1=ztrans(f1,n,z); pretty(F1)

inv_F1=iztrans(F1,z,n)

z ----- z - 1 inv_F1 = 1

%

<11>

（4）

clear

syms n z clear % delta=sym('kroneckerDelta(n, 0)'); % KD=ztrans(delta,n,z) inv_KD=iztrans(KD) KD = 1

inv_KD =

kroneckerDelta(n, 0)

<16> <17>

（5）

syms n z clear % <20> k=sym('k','positive'); % <21> fd=sym('f(n)*kroneckerDelta(n-k, 0)'); FD=ztrans(fd,n,z)

inv_FD=iztrans(FD,z,n) FD =

piecewise([k in Z_, f(k)/z^k], [Otherwise, 0]) inv_FD =

piecewise([k in Z_, f(k)*kroneckerDelta(k - n, 0)], [Otherwise, 0])

（6）

syms a z n clear GZ=exp(-a/z); gn=iztrans(GZ,z,n) gn =

(-a)^n/factorial(n)

% %

<28>

25

2.5.4

【例2.5-6】

符号卷积

syms T t tao

ut=exp(-t); % ht=exp(-t/T)/T; % uh_tao=subs(ut,t,tao)*subs(ht,t,t-tao); % yt=simple(simple(int(uh_tao,tao,0,t))) % yt =

-(1/exp(t) - 1/exp(t/T))/(T - 1)

【例2.5-7】

syms s

yt=ilaplace(laplace(ut,t,s)*laplace(ht,t,s),s,t); yt=simple(yt) yt =

-(1/exp(t) - 1/exp(t/T))/(T - 1)

2.6

2.6.1

符号矩阵分析和代数方程解

符号矩阵分析

【例2.6-1】 （1）

syms a11 a12 a21 a22 A=[a11,a12;a21,a22] DA=det(A) IA=inv(A) A =

[ a11, a12] [ a21, a22] DA =

a11*a22 - a12*a21 IA =

[ a22/(a11*a22 - a12*a21), -a12/(a11*a22 - a12*a21)] [ -a21/(a11*a22 - a12*a21), a11/(a11*a22 - a12*a21)]

（2）

EA=subexpr(eig(A),'D') D =

(a11^2 - 2*a11*a22 + a22^2 + 4*a12*a21)^(1/2) EA =

a11/2 + a22/2 - D/2 a11/2 + a22/2 + D/2

【例2.6-2】 （1）

syms t

A=sym([sqrt(3)/2,1/2;1/2,sqrt(3)/2]) G=[cos(t),-sin(t);sin(t),cos(t)]; GA=G*A A =

[ 3^(1/2)/2, 1/2] [ 1/2, 3^(1/2)/2]

%

26

GA =

[ (3^(1/2)*cos(t))/2 - sin(t)/2, cos(t)/2 - (3^(1/2)*sin(t))/2] [ cos(t)/2 + (3^(1/2)*sin(t))/2, sin(t)/2 + (3^(1/2)*cos(t))/2]

（2）

clf

An=subs(GA,t,110/180*pi); % Op=[0;0]; % Ad=double(A); % v1=[Op,Ad(:,1)]';v2=[Op,Ad(:,2)]'; % u1=[Op,An(:,1)]';u2=[Op,An(:,2)]'; %

plot(v1(:,1),v1(:,2),'--k',v2(:,1),v2(:,2),'b') axis([-1,1,-1,1]),axis square hold on

hu=plot(u1(:,1),u1(:,2),'--k',u2(:,1),u2(:,2),'b'); set(hu,'LineWidth',4) title('Givens Rotation')

Lstr=['旋转前的 v1';'旋转前的 v2';'旋转后的 u1';'旋转后的 u2']; % <16> legend(Lstr,'Location','South') % <17> hold off grid on Givens Rotation10.80.60.40.20-0.2-0.4-0.6-0.8-1 -1-0.5旋转前的 v1旋转前的 v2旋转后的 u1旋转后的 u200.51 图 5.6-1 Givens旋转的几何意义

2.6.2

【例 2.6-3】 （1）

线性方程组的符号解

A=sym([1 1/2 1/2 -1;1 1 -1 1;1 -1/4 -1 1;-8 -1 1 1]); b=sym([0;10;0;1]); X1=A\\b % X1 = 1 8 8 9

27

（2）

eq1=sym('d+n/2+p/2-q'); % <4> eq2=sym('d+n-p+q-10'); % eq3=sym('d-n/4-p+q'); % eq4=sym('-8*d-n+p+q-1'); %

S=solve(eq1,eq2,eq3,eq4,'d','n','p','q'); % <8> disp([' d',' n',' p',' q']) disp([S.d,S.n,S.p,S.q]) % <10> d n p q [ 1, 8, 8, 9]

2.6.3

一般代数方程组的解

【例2.6-4】

S=solve('u*y^2+v*z+w=0','y+z+w=0','y','z') % <1>

disp('S.y'),disp(S.y),disp('S.z'),disp(S.z) % <2> S =

y: [2x1 sym] z: [2x1 sym] S.y

(v + 2*u*w + (v^2 + 4*u*w*v - 4*u*w)^(1/2))/(2*u) - w (v + 2*u*w - (v^2 + 4*u*w*v - 4*u*w)^(1/2))/(2*u) - w S.z

-(v + 2*u*w + (v^2 + 4*u*w*v - 4*u*w)^(1/2))/(2*u) -(v + 2*u*w - (v^2 + 4*u*w*v - 4*u*w)^(1/2))/(2*u)

【例2.6-5】

syms d n p q

eq1=d+n/2+p/2-q;eq2=n+d+q-p-10;eq3=q+d-n/4-p; S=solve(eq1,eq2,eq3,d,n,p,q);

disp([' S.d',' S.n',' S.p',' S.q']) disp([S.d,S.n,S.p,S.q])

S.d S.n S.p S.q [ z/3 - 2, 8, (4*z)/3 - 4, z]

【例2.6-6】

clear all,syms x;

s=solve('(x+2)^x=2','x') xs=(s(1)+2)^s(1) % 验算

s =

matrix([[0.69829942170241042826920133106081]]) xs = 2.0

2.7

代数状态方程求符号传递函数

2.7.1

结构框图的代数状态方程解法

【例2.7-1】

28

% <2> % <3>

<3>

（1） （2）

（3）

（4）

syms G1 G2 G3 G4 H1 H2 H3

A=[ 0, 0, 0, 0, 0, 0, -G1; G2, 0, 0, 0, 0, -G2, 0; 0, G3, 0, 0, G3, 0, 0; 0, 0, G4, 0, 0, 0, 0; 0, 0, 0, H2, 0, 0, 0; 0, 0, 0, H1, 0, 0, 0; 0, 0, 0, H3, 0, 0, 0]; b=[ G1; 0; 0; 0; 0; 0; 0]; c=[ 0, 0, 0, 1, 0, 0, 0];

Y2Ua=c*((eye(size(A))-A)\\b); %利用“左除”取代“求逆”，计算传递函数 disp([blanks(5),'传递函数 Y2Ua 为']) pretty(Y2Ua)

传递函数 Y2Ua 为

G1 G2 G3 G4

------------------------------------------- G2 G3 G4 H1 - G3 G4 H2 + G1 G2 G3 G4 H3 + 1

（2.7-2） （2.7-4）

2.7.2

【例2.7-2】 （1）

信号流图的代数状态方程解法

此

（2） （3）

syms G1 G2 G3 G4 H1 H2 H3 A=[ 0, 0, 0, 0, -H3; G1, 0, 0, 0, -H1; 0, G2, 0, 0, H2; 0, 0, G3, 0, 0; 0, 0, 0, G4, 0]; b=[ 1; 0; 0; 0; 0]; c=[ 0, 0, 0, 0, 1]; Y2Ub=c*((eye(size(A))-A)\\b); %

disp([blanks(5),'传递函数 Y2Ub 为']) pretty(Y2Ub)

传递函数 Y2Ub 为

G1 G2 G3 G4

------------------------------------------- G2 G3 G4 H1 - G3 G4 H2 + G1 G2 G3 G4 H3 + 1

（2.7-5） （2.7-6）

（4）

syms s

Sblock={100/(s+10),1/(s+1),(s+1)/(s^2+4*s+4),(s+1)/(s+6),(2*s+12)/(s+1),(s+1)/(s+2),1}; % <9> ww=subs(Y2Ub,{G1,G2,G3,G4,H1,H2,H3},Sblock);% <10> Y2Uc=simple(ww);

[NN,DD]=numden(Y2Uc); %

29

NN=expand(NN); % disp('参数具体化的传递函数 Y2Uc 为') pretty(NN/DD)

参数具体化的传递函数 Y2Uc 为

2

100 s + 300 s + 200

------------------------------------------- 5 4 3 2

s + 21 s + 157 s + 663 s + 1301 s + 910

?

2.8

符号计算结果的可视化

直接可视化符号表达式 单独立变量符号函数的可视化

2.8.1 1

【例2.8-1】

syms t tao

y=2/3*exp(-t/2)*cos(sqrt(3)/2*t) % s=subs(int(y,t,0,tao),tao,t) % subplot(2,1,1) % ezplot(y,[0,4*pi]),ylim([-0.2,0.7]) grid on % subplot(2,1,2) % ezplot(s,[0,4*pi]) grid on

title('s = \\int y(t)dt') % y =

(2*cos((3^(1/2)*t)/2))/(3*exp(t/2)) s =

1/3-(2*cos((3^1/2)*t)/2)/2-3^(1/2)*sin((3^(1/2)*t)/2))/2))/(3*exp(t/2))

30

(2 cos((31/2 t)/2))/(3 exp(t/2))0.60.40.20-0.20246ts = ? y(t)dt0.50.40.30.20246t8101281012图 2.8-1 ezplot使用示例

2 双独立变量符号函数的可视化

【例2.8-2】

clf

x='cos(s)*cos(t)'; % y='cos(s)*sin(t)'; % z='sin(s)'; % ezsurf(x,y,z,[0,pi/2,0,3*pi/2]) % view(17,40) % shading interp colormap(spring) %

light('position',[0,0,-10],'style','local') light('position',[-1,-0.5,2],'style','local') material([0.5,0.5,0.5,10,0.3])

% %

31

图 2.8-2 ezsurf在参变量格式下绘制的图形

2.8.2

【例2.8-3】 （1）

符号计算结果的数值化绘图

clear

syms x y real % fx=1-2/(1+exp(x)); % disp('f(x)=') pretty(fx) disp(' ')

fxint=int(fx,x,0,x) % f(x)=

2

1 - ---------- exp(x) + 1

fxint =

log((exp(x) + 1)^2/4) - x

（2）

xk=0:0.1:2; % fxk=subs(fx,x,xk); % fxintk=subs(fxint,x,xk); %

plot(xk,fxk,'g',xk,fxintk,'r','LineWidth',2.5) % title

xlabel('x')

legend('f(x)','\\int^x_0 f(x) dx','Location','best')

32

函数及其积分函数0.90.80.70.60.50.40.30.20.10 0f(x) ?x f(x) dx00.20.40.60.81x1.21.41.61.82图 2.8-3 用双精度数据绘曲线

（3）

gy=subs(finverse(fx),x,y) % gyint=int(gy,y,0,y) % gy =

log(-(y + 1)/(y - 1)) gyint =

piecewise([y < 1, log(1 - y^2) + y*log(y + 1) - y*log(1 - y)], [1 <= y, log(y^2 - 1) + y*log(-(y + 1)/(y - 1)) - pi*i])

（4）

gf=simplify(subs(gy,y,fx)) gf = x

%

（5）

yk=subs(fx,x,xk); % gyintk=subs(gyint,y,yk); % GYintk=xk.*fxk-fxintk; %

plot(yk,gyintk,'r','LineWidth',5) hold on

plot(yk,GYintk,'+k') hold off xlabel('y')

legend('直接法计算反函数积分', '互补法求反函数积分','location','best')

33

0.7直接法计算反函数积分互补法求反函数积分 0.60.50.40.30.20.10 00.10.20.30.4y0.50.60.70.8图 2.8-4 反函数积分两种算法结果比较

2.8.3

【例2.8-4】 （1）

可视化与数据探索

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

（2）

Fxy=sym('sin(x^2+y)')

Fxy_TL1=Fxy-TL1 % Fxy =

sin(x^2 + y) Fxy_TL1 =

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

（3）

figure(1) % ezsurf(Fxy,[-2,2,-3,3]) % shading interp % view([-63,52]) % colormap(spring)

light,light('position',[-10,4,50],'style','local','color','r')

34

图 2.8-5 原函数在较大范围内的图形

（4）

figure(2)

ezsurf(TL1,[-2,2,-3,3]) shading interp view([-43,54]) colormap(spring) light

light('position',[-10,2,2],'style','local','color',[0.8,0.3,0.3]) light('position',[-2,-10,2],'style','local','color',[0.4,0.5,0.7])

35

图 2.8-6 Taylor展开在较大范围内的图形

（5）

figure(3)

ezsurf(Fxy,[-0.5,0.5,-0.5,0.5],'circ') % axis([-1,1,-1,1,-2,2]) shading interp colormap(spring) view([-49,17]) light

light('position',[-30,0,-2],'style','local','color','r')

2'r'

图 2.8-7 原函数在小范围内的图形

（6）

figure(4)

ezsurf(Fxy_TL1,[-0.5,0.5],'circ') shading interp colormap(spring) view([-53,34]) light

light('position',[-10,15,0],'style','local','color',[0.8,0.3,0.3])

36

图 2.8-8 小范围内的误差曲面

2.9

符号计算资源深入利用

符号表达式、串操作及数值计算M码间的转换 符号工具包资源表达式转换成M码函数 把符号包资源转换成M码函数的示例

2.9.1 2.9.2

1

【例2.9-1】 （1）

s=dsolve('x*D2y-3*Dy=x^2','y(a)=yL,y(b)=yR','x') s =

(a^4*b^3 + 3*yR*a^4 - a^3*b^4 - 3*yL*b^4)/(3*a^4 - 3*b^4) - x^3/3 + (x^4*(a^3 - b^3 + 3*yL - 3*yR))/(3*a^4 - 3*b^4)

（2）

Hs=matlabFunction(s,'file','exm050901_ZZY','vars',{'x','a','b','yL','yR'},'outputs',{'y'}) % <2> Hs =

@exm050901_ZZY

（3）

a=1;b=5;yL=0;yR=0;

xn=-1:6;yn=Hs(xn,a,b,yL,yR) x=-1:0.2:6;y=Hs(x,a,b,yL,yR); plot(x,y,'b-'),hold on

% %

<4> <5>

37

plot([1,5],[0,0],'.r','MarkerSize',20),hold off

title(['xy{\\prime\\prime}- 3y{\\prime}=x^2,',' ','y(1)=0, y(5)=0']) text(1,1,'y(1)=0'),text(4,1,'y(5)=0') xlabel('x'),ylabel('y') yn =

Columns 1 through 7

0.6667 0.2671 0.0000 -1.3397 -3.3675 -4.1090 0.0000 Column 8 14.1132 xy??- 3y?=x2, y(1)=0, y(5)=015105yy(1)=00y(5)=0-5-1012x3456图5.9-2 两点边值问题的解曲线

（4）

function y = exm050901_ZZY(x,a,b,yL,yR) % % % % % % % % %

t7 = a.^2; t8 = b.^2; t9 = t7.^2; t10 = t8.^2; t11 = x.^2; t12 = 3.*t10; t13 = t12 - 3.*t9; t14 = 1./t13;

y = t14.*(3.*t10.*yL - 3.*t9.*yR + a.*t10.*t7 - b.*t8.*t9) - (x.*t11)./3 - t14.*(3.*yL - 3.*yR +

38

a.*t7 - b.*t8).*t11.^2;

（5）

ym=exm050901_ZZY(xn,a,b,yL,yR) % <12> ym =

Columns 1 through 7

0.6667 0.2671 0.0000 -1.3397 -3.3675 -4.1090 0.0000 Column 8 14.1132

? 〖

【例2.9-2】

（1）

function y = exm050902_ZZY(de,x,a,b,yL,yR,flag) % % % % % % % %

if nargin~=7

error('输入量数目应有 7 个 !') end

s=dsolve(de,'y(a)=yL,y(b)=yR','x');

Hs=matlabFunction(s,'vars',{'x','a','b','yL','yR'}); y=Hs(x,a,b,yL,yR); if flag==1

plot(x,y,'-b',[a,b],[yL,yR],'*r') title(de)

xlabel('x'),ylabel('y') shg end

（2）

x=-1:6;

a=1;b=5;yL=0;yR=0; de1='x*D2y-3*Dy=x^2';

y = exm050902_ZZY(de1,x,a,b,yL,yR,0) y =

0.6667 0.2671 0.0000 -1.3397 -3.3675 -4.1090 0 14.1132

（3）

de2='x*D2y-3*Dy=3*x^2+x'; x=-1:0.2:6;

y = exm050902_ZZY(de2,x,3,6,-2,10,1);

39

% <14>

x*D2y-3*Dy=3*x2+x151050y-5-10-15-1012x3456图5.9-3 用exm050902_ZZY求解另一微分方程所得解

2.9.3

【例2.9-3】 （1）

借助mfun调用MuPAD特殊函数

syms t

gt=1/log(t);

gt_0=limit(gt, t,0,'right') % gt_0 = 0

（2）

ezplot(gt,[0,1]) grid on

legend('gt')

%

40

1/log(t) 0-2-4-6-8-10-12-14gt 00.10.20.30.40.5t0.60.70.80.91图5.10-1 在关心区间内的被积函数

（3）

fx=int(gt,t,0,x) % <7> Warning: Explicit integral could not be found. fx =

piecewise([x < 1, Li(x)], [Otherwise, int(1/log(t), t = 0..x)])

（4）

x=0:0.05:0.9; % <8>

fxMfun=mfun('Li',x) % <9> % fxMfun =

Columns 1 through 8

NaN -0.0131 -0.0324 -0.0564 -0.0851 -0.1187 -0.1574 -0.2019

Columns 9 through 16

-0.2529 -0.3114 -0.3787 -0.4564 -0.5469 -0.6534 -0.7809 -0.9369

Columns 17 through 19

-1.1340 -1.3959 -1.7758

（5）

hold on

plot(x,fxMfun,'--r','LineWidth',3)

legend('gt','fxMfun','Location','Best') title(' ') % hold off

41

0-2-4-6-8-10-12-14 0gtfxMfun0.10.20.30.40.5t0.60.70.80.91图5.10-2 被积函数曲线gt和积分曲线fx

（6）

fx_matlab=-expint(-log(x)) % <15> % fx_matlab =

Columns 1 through 8

NaN -0.0131 -0.0324 -0.0564 -0.0851 -0.1187 -0.1574 -0.2019

Columns 9 through 16

-0.2529 -0.3114 -0.3787 -0.4564 -0.5469 -0.6534 -0.7809 -0.9369

Columns 17 through 19

-1.1340 -1.3959 -1.7758

习题2

1． 说出以下四条指令产生的结果各属于哪种数据类型，是“双精度”对象，还是“符号”

对象？

3/7+0.1, sym(3/7+0.1), vpa(sym(3/7+0.1))

2． 在不加专门指定的情况下，以下符号表达式中的哪一个变量被认为是独立自由变量。

sym('sin(w*t)') , sym('a*exp(-X)' ) , sym('z*exp(j*th)') 3． 求以下两个方程的解： （提示：关于符号变量的假设要注意）

（1）试写出求三阶方程x?44.5?0正实根的程序。注意：只要正实根，不要出现其

他根。

（2）试求二阶方程x?ax?a?0在a?0时的根。

4． 观察一个数（在此用@记述）在以下四条不同指令作用下的异同：

a = @ , b = sym( @ ), c = sym( @ ,'d ' ), d = sym( '@ ' )

在此，@ 分别代表具体数值 7/3 , pi/3 , pi*3^(1/3) ；而异同通过vpa(abs(a-d)) , vpa(abs(b-d)) , vpa(abs(c-d))等来观察。

223 42

?a11?5． 求符号矩阵A?a21???a31简洁化。

a12a22a32a13?a23??的行列式值和逆，所得结果应采用“子表达式置换”a33???kf(k)z???a6． 对函数f(k)???0?kk?0，当a为正实数时，求k?0k?02k?1。（提示：symsum。

实际上，这就是根据定义求Z变换问题。）

2?x?1?7． 对于x?0，求???2k?1x?1??k?0。（提示：理论结果为lnx；注意限定性假设）

dydy8． （1）通过符号计算求y(t)?sint的导数。（2）然后根据此结果，求

dtdtdydtt??9． 求出

。

和

t?0??21.7??5?e?x2sinxdx的具有64位有效数字的积分值。（提示：int, vpa, ezplot）

x2110． 计算二重积分

??1(x2?y2)dydx。

11． 在[0,2?]区间，画出y(x)?12． 在n?0的限制下，求

?x0sintdt曲线，并计算y(4.5)。（提示：int, subs） t?20y(n)??sinnxdx的一般积分表达式，并计算y()的

1332位有效数字表达。（提示：注意限定条件；注意题目要求32位有效）

kk13． 有序列x(k)?a，h(k)?b，（在此k?0，a?b），求这两个序列的卷积

y(k)??h(n)x(k?n)。（提示：symsum, subs）

n?0k?3t14． 设系统的冲激响应为h(t)?e，求该系统在输入u(t)?cost，t?0作用下的输出。

（提示：直接卷积法，变换法均可；）

15． 求f(t)?Ae??t,??0的Fourier变换。（提示：注意限定）

????t??的Fourier变换，并画出A?2,??2时的幅频谱。（提t????t?A?1?16． 求f(t)??????0?17． 求F(s)?示：注意限定；simple）

s?3的Laplace反变换。

s3?3s2?6s?418． 利用符号运算证明Laplace变换的时域求导性质：L?示：用sym('f(t)')定义函数f(t)） 19． 求

?df(t)?（提??s?L?f(t)??f(0)。dt??f(k)?ke?? k T的Z变换表达式。

2220． 求方程x?y?1,xy?2的解。（提示：正确使用solve）

21． 求图p2-1所示信号流图的系统传递函数，并对照胡寿松主编“自动控制原理”中的例

43

2-21结果，进行局部性验证。（提示：在局部性验证时，把不存在支路的符号增益设置为0）

图p2-1

22． 采用代数状态方程法求图p2-2所示结构框图的传递函数

YY和。（提示：列出正确UW?x?Ax?bU?fW的状态方程?，进而写出相关输入输出对之间的传递函数表达式

Y?cx?dU?gW?YY?c(I?A)?1b?d和?c(I?A)?1f?g） UN图p2-2

23． 求微分方程

yy??x?0的通解，并绘制任意常数为1时，如图p2-3所示的解曲线

54图形。（提示：通解中任意常数的替代；构造能完整反映所有解的统一表达式，然后

绘图。）

44

图 p2-3 微分方程的解曲线

24． 求一阶微分方程x??at2?bt,x(0)?2的解。 25． 求边值问题

dfdg?3f?4g,??4f?3g,f(0)?0,g(0)?1的解。 dxdx 45

本文来源：https://www.bwwdw.com/article/t8w8.html