ch2 - 符号计算2010a(11)

图 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