ch2 - 符号计算2010a(5)

（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