%习题2-1
class（vpa（sym（3/7+0.1）,4））
%习题2-2
findsym（sym（'sin（w*t）'）,1）
findsym（sym（'a*exp（-X）'）,1）
findsym（sym（'z*exp（j*th）'）,1）
%习题2-3
symsxapositive
x_1=solve（'x^3-44.5'）
symsxunreal
x_2=solve（'x^2-a*x+a^2'）
symsaunreal
x_2=solve（'x^2-a*x+a^2'）
x_3=solve（'x^2-a*x-a^2'）
symsxreal
evalin（symengine,'anames（Properties）'）
evalin（symengine,'getprop（x）'）
%习题2-4
7/3,pi/3,pi*3^（1/3）
a=pi*3^（1/3）,
b=sym（pi*3^（1/3））,
c=sym（pi*3^（1/3）,'d'）,
d=sym（'pi*3^（1/3）'）
vpa（abs（a-d））,vpa（abs（b-d））,vpa（abs（c-d））
%习题2-5
symsa11a12a13a21a22a23a31a32a33
A=[a11a12a13;a21a22a23;a31a32a33]
a=det（A）
B=inv（A）
C=subexpr（B）
[RS,w]=subexpr（B,'w'）
%习题2-6
symsk
symsapositive
%fk=a^k*heaviside（k）
fk=a^k
s=symsum（fk,k,0,inf）
%习题2-7
clearall
symsk
symsxpositive
fk=2/（2*k+1）*（（（x-1）/（x+1））^（2*k+1））
s=symsum（fk,k,0,inf）
s1=simple（s）
%习题2-8
clearall,symst
y=abs（sin（t））
df=diff（y）,class（df）
df1=limit（df,t,0,'left'）
df2=subs（df,'t',sym（pi/2））
%习题2-9
clearall,symsx;
f=exp（-abs（x））.*abs（sin（x））;
fint=int（f,x,-5*pi,1.7*pi）,digits（64）,vpa（fint）
class（fint）
%习题2-10
clearall,symsxy,f=x.^2+y.^2,
fint=（int（int（f,y,1,x.^2）,x,1,2））,double（fint）
%习题2-11
clearall,symstx;f=sin（t）/t,yx=int（f,t,0,x）,ezplot（yx,[02*pi]）
fint=subs（yx,x,4.5）,%或yxd=int（f,t,0,4.5）,fint=double（yxd）
holdon,plot（4.5,fint,'*r'）
%习题2-12
%clearall,symsxn;f=（sin（x））^n;yn=int（f,'x',0,pi/2）,class（yn）
clearall,symsxn;symsnpositive;f=（sin（x））^n;yn=int（f,'x',0,pi/2）,class（yn）
%y（1/3）=?
yn1=subs（yn,'n',sym（1/3））,vpa（yn1）
%或yn=limit（yn,n,1/3）,vpa（yn）
%或yy=int（sin（x）.^（1/3）,x,0,pi/2）,vpa（yy）
%习题2-23
clear,symsxyS
S=dsolve（'Dy*y/5+x/4=0','x'）
ezplot（subs（S
（1）,'C3',1）,[-2,2-2,2],1）,holdon
ezplot（subs（S
（2）,'C3',1）,[-2,2-2,2],1）
%解为S=
%2^（1/2）*（C3-（5*x^2）/8）^（1/2）
%-2^（1/2）*（C3-（5*x^2）/8）^（1/2）
ezplot（subs（S
（1）,'C3',1）,[-2,2-2,2],1）,holdon
ezplot（subs（S
（2）,'C3',1）,[-2,2-2,2],1）%用此两条指令绘圆，在y=0处有间隙
ezplot（subs（y^2-（S
（1））^2,'C3',1）,[-2,2-2,2],2）%用椭圆方程绘图不产生间隙
colormap（[001]）%用ezplot（fun）绘图时，如果fun中只有一个参数，绘图的颜色是蓝色；如果fun中有两个参数，绘图的颜色是绿色，此指令设置图形颜色为蓝。
gridon
S1=subs（S
（1）,'C3',1）
subs（S1,{x},{1.6^（1/2）}）
y=double（solve（S1））
t=linspace（y
（2）,y
（1）,100）
S2=subs（S
（2）,'C3',1）
figure
plot（t,subs（S1,x,t））,holdon
plot（t,subs（S2,x,t））
axis（[-2,2-2,2]）
%习题2-24
x=dsolve（'Dx=a*t^2+b*t','x（0）=2','t'）
%习题2-25
[f,g]=dsolve（'Df=3*f+4*g','Dg=-4*f+3*g','f（0）=0','g（0）=1','x'）
第三章
%习题3-1
a=0:
2*pi/9:
2*pi
b=linspace（0,2*pi,10）
%习题3-2
rand（'twister',0）,
A=rand（3,5）
[I1,J1]=find（A>0.5）
subindex_A=sub2ind（size（A）,I1,J1）
subindex_A=find（A>0.5）
[I,J]=ind2sub（size（A）,subindex_A）
index_A=[IJ]
%122_2
rand（'twister',0）,A=rand（3,5）,B=A>0.5,C=find（B）
[ci,cj]=ind2sub（size（A）,C）%此法太繁
[Ci,Cj]=find（B）
%122_5
t=linspace（1,10,100）
（1）y=1-exp（-0.5.*t）.*cos（2.*t）,plot（t,y）
（2）L=length（t）
fork=1:
L,yy（k）=1-exp（-0.5*t（k））*cos（2*t（k））,end,plot（t,yy）
%122_6
clear,formatlong,rand（'twister',1）,A=rand（3,3）,
B=diag（diag（A））,C=A.*（~B）%或C=A-B%或C=triu（A,1）+tril（A,-1）
%习题3-3
%s=sign（randint（1,1000,[],123）-.5）;
%n=sum（s==-1）
rand（'state',123）,
s=sign（rand（1,1000）-.5）,
n=sum（s==-1）
%习题3-4
clearall
A=[12;34]
B1=A.^（0.5）,B2=A^（0.5）,
A1=B1.^2
A2=B2^2
A1-A
norm（A1-A）
A1-A2
norm（A1-A2）
clearall
A=sym（'[12;34]'）
B1=A.^（0.5）,B2=A^（0.5）,
A1=simple（B1.^2）
A2=simple（B2^2）
A1-A
vpa（A1-A）
A1-A2
vpa（A1-A2）
%习题3-5
t=linspace（1,10,100）
%
（1）
L=length（t）
fork=1:
L
yy（k）=1-exp（-0.5*t（k））*cos（2*t（k））
end
figure
（1）,plot（t,yy）
%
（2）
y=1-exp（-0.5.*t）.*cos（2.*t）,
figure
（2）,plot（t,y）
figure（3）,subplot（1,2,1）,plot（t,y）
subplot（1,2,2）,plot（t,yy）
%习题3-6
clear,formatlong,
rand（'twister',1）,A=rand（3,3）,
B=diag（diag（A））,
C=A.*（~B）
%或C=A-B%或C=triu（A,1）+tril（A,-1）
%习题3-7
clear,
x=-3*pi:
pi/15:
3*pi;y=x;
[X,Y]=meshgrid（x,y）;%X=ones（size（y））*x;Y=y*ones（size（x））;
warningoff;
Z=sin（X）.*sin（Y）./X./Y;
%
（1）“非数”数目
m=sum（sum（isnan（Z）））;%k=Z（isnan（Z））;m=length（k）
%
（2）绘图
surf（X,Y,Z）;shadinginterp
%（3）“无裂缝”图形的全部指令：
x=-3*pi:
pi/15:
3*pi;y=x';
[X,Y]=meshgrid（x,y）;%X=ones（size（y））*x;Y=y*ones（size（x））;
X=X+（X==0）*eps;
Y=Y+（Y==0）*eps;
Z=sin（X）.*sin（Y）./X./Y;
surf（X,Y,Z）;shadinginterp
%习题3-8
clear
fork=10:
-1:
1;
A=reshape（1:
10*k,k,10）;
Sa（k,:
）=sum（A,1）;
end
Sa
clear
fork=10:
-1:
1;
A=reshape（1:
10*k,k,10）;
ifk==1
Sa（k,:
）=A;
else
Sa（k,:
）=sum（A）;
end
end
Sa
第四章
%习题4-1
clearall
loadprob_data401;
diff_y=diff（y）./diff（t）;
gradient_y=gradient（y）./gradient（t）;
%plot（t（1:
end-1）,diff_y,t,gradient_y）
figure
（1）
subplot（1,2,1）,plot（t,y,t（1:
end-1）,diff_y）
subplot（1,2,2）,plot（t,y,t,gradient_y）
%上面结果不好因自变量采样间距太小，将间距增大后结果较好
N=20
diff_y1=（diff（y（1:
N:
end）））./diff（t（1:
N:
end））;
gradient_y1=（gradient（y（1:
N:
end）））./gradient（t（1:
N:
end））;
%plot（t（1:
end-1）,diff_y,t,gradient_y）
t1=t（1:
N:
end）;
length（t1）
figure
（2）
subplot（1,2,1）,plot（t,y,t1（1:
end-1）,diff_y1）
subplot（1,2,2）,plot（t,y,t1,gradient_y1）
%习题4-2
d=0.5;tt=0:
d:
10;t=tt+（tt==0）*eps;y=sin（t）./t;s=d*trapz（y）%计算出积分值
ss=d*（cumtrapz（y））%计算梯形法累计积分并绘积分曲线
plot（t,y,t,ss,'r'）,holdon
%用find指令计算y（4.5）,并绘出该点
y4_5=ss（find（t==4.5））
%插值法计算y（4.5）,并绘出该点
yi=interp1（t,ss,4.5）,plot（4.5,yi,'r+'）
%yi=interp1（t,ss,[4.24.34.5]）,plot（[4.24.34.5],yi,'r+'）
%clf
%用精度可控指令quad即Simpson法,quadl即Lobatto法计算y（4.5）
yy=quad（'sin（t）./t',0,4.5）
yy=quadl（'sin（t）./t',0,4.5）
warningoff%取消警告性提示时用
%此法可用，但有警告性提示,取消提示加warningoff
clear
tt=0:
0.1:
10;warningoff
fori=1:
101
q（i）=quad（'sin（t）./t',0,tt（i））;
end
plot（tt,q）
y=quad（'sin（t）./t',0,4.5）
%符号解法，匿名函数求y（4.5）
f=@（x）（int（'sin（t）/t',0,x））,vpa（f（4.5））
%符号解法
symsxty1y2y1i,y1=sin（t）./t,y1i=int（y1,t,0,x）,y2=subs（y1i,x,4.5）
holdon,plot（4.5,y2,'*m'）
%习题4-3
clearall
d=pi/20;x=0:
d:
pi;fx=exp（sin（x）.^3）;
s=d*trapz（fx）%梯形法计算积分值
%用精度可控指令quad即Simpson法,quadl即Lobatto法计算积分
s1=quad（'exp（sin（x）.^3）',0,pi）
s2=quadl（'exp（sin（x）.^3）',0,pi）
%符号计算解法
%s3=vpa（int（'exp（sin（x）^3）','x',0,pi））
s3=vpa（int（'exp（sin（x）^3）',0,pi））
s4=vpa（int（sym（'exp（sin（x）^3）'）,0,pi））
%习题4-4
clearall
%用精度可控指令quad即Simpson法,quadl即Lobatto法计算积分
s1=quad（'exp（-abs（x））.*abs（sin（x））',-5*pi,1.7*pi,1e-10）
s2=quadl（'exp（-abs（x））.*abs（sin（x））',-5*pi,1.7*pi）
%符号计算解法
在Matlab6.5版本上下式计算结果多加了值1
%s3=vpa（int（'exp（-abs（x））*abs（sin（x））',-5*pi,1.7*pi））
symsx;
s3=vpa（int（exp（-abs（x））*abs（sin（x））,-5*pi,1.7*pi））
%s3=vpa（int（'sin（x）',-pi/2,pi/2））
%梯形法计算积分值
d=pi/1000;x=-5*pi:
d:
1.7*pi;fx=exp（-abs（x））.*abs（sin（x））;
s=d*trapz（fx）
%习题4-5
clearall
x1=-5;x2=5;
%采用内联函数或字符串函数求极值
yx=inline（'（sin（5*t））.^2.*exp（0.06*t.^2）-1.5.*t.*cos（2*t）+1.8.*abs（t+0.5）'）
[xn0,fval]=fminbnd（yx,x1,x2）
%绘函数图并标出最小值点
t=x1:
0.1:
x2;plot（t,yx（t））,holdon,plot（xn0,fval,'r*'）
%字符串函数求极值
[xn0,fval]=fminbnd（'（sin（5.*x））.^2.*exp（0.06.*x.^2）-1.5.*x.*cos（2.*x）+1.8.*abs（x+0.5）',-5,5）
%yx='（sin（5.*x））.^2.*exp（0.06.*x.^2）-1.5.*x.*cos（2.*x）+1.8.*abs（x+0.5）'
%[xn0,fval]=fminbnd（yx,x1,x2）
下一条指令即字符串函数在无论在2010a版本还是Matlab6.5版本上不适用，因为对字符串函数只是别符号x,不识别t
%[xn0,fval]=fminbnd（'（sin（5.*t））.^2.*exp（0.06.*t.^2）-1.5.*t.*cos（2.*t）+1.8.*abs（t+0.5）',-5,5）
下一条指令即匿名函数在Matlab6.5版本上不适用
%yx=@（t）（（sin（5*t））.^2.*exp（0.06*t.^2）-1.5.*t.*cos（2*t）+1.8.*abs（t+0.5））%本条指令即匿名函数在Matlab6.5版本上不适用
%yx=@（t）（（sin（5.*x））.^2.*exp（0.06.*x.^2）-1.5.*x.*cos（2.*x）+1.8.*abs（x+0.5））%本条指令即匿名函数在Matlab6.5版本上不适用
%[xn0,fval]=fminbnd（yx,x1,x2）
%习题4-6
tspan=[0,0.5];%
y0=[1;0];%
[tt,yy]=ode45（@DyDt_6,tspan,y0）;
y0_5=yy（end,1）
%figure
（1）
%plot（tt,yy（:
1））
%xlabel（'t'）,title（'x（t）'）
%figure
（2）
%plot（yy（:
1）,yy（:
2））%
%xlabel（'位移'）,ylabel（'速度'）
函数DyDt_6另存为一个文件
functionydot=DyDt_6（t,y）
mu=3;
ydot=[y
（2）;mu*y
（2）-2*y
（1）+1];
%符号计算解法
S=dsolve（'D2y-3*Dy+2*y=1','y（0）=1','Dy（0）=0'）
ys0_5=subs（S,0.5）
%习题4-7
clearall
A=magic（8）
B=orth（A）
rref（A）
rref（B）
A=
642361606757
955541213515016
1747462021434224
4026273736303133
3234352928383925
4123224445191848
4915145253111056
858595462631
B=
-0.353553390593270.540061724867320.35355339059327
-0.35355339059327-0.38575837490523-0.35355339059327
-0.35355339059327-0.23145502494314-0.35355339059327
-0.353553390593270.077151674981050.35355339059327
-0.35355339059327-0.077151674981050.35355339059327
-0.353553390593270.23145502494314-0.35355339059327
-0.353553390593270.38575837490523-0.35355339059327
-0.35355339059327-0.540061724867320.35355339059327
ans=
10011001
01034-3-47
001-3-445-7
00000000
00000000
00000000
00000000
00000000
%习题4-8
A=sym（gallery（5））
[v,lambda]=eig（A）
cond（A）
clearall,
A=gallery（5）
[V,D,s]=condeig（A）
[V,D]=eig（A）
cond（A）
jordan（A）
ans=
01000
00100
00010
00001
00000
%习题4-9
clearall,
A=magic（3）
b=ones（3,1）
x=A\b
x=
0.06666666666667
0.06666666666667
0.06666666666667
x=inv（A）*b
rref（[A,b]）
ans=
1.00000000000000000.06666666666667
01.0000000000000000.06666666666667
001.000000000000000.06666666666667
%习题4-10
解不唯一
clearall,
A=magic（4）
b=ones（4,1）
rref（[A,b]）
ans=
1.00000000000000001.000000000000000.05882352941176
01.0000000000000003.000000000000000.11764705882353
001.00000000000000-3.00000000000000-0.05882352941176
00000
x=A\b
x=
0.05882352941176
0.11764705882353
-0.05882352941176
0
xg=null（A）
xg=
0.22360679774998
0.67082039324994
-0.67082039324994
-0.22360679774998
x+xg也是方程的解
ans=
0.28243032716174
0.78846745207347
-0.72964392266170
-0.22360679774998
升级会员 