z2=a*ones(size(x));
subplot(1,3,2),mesh(x,y,z2),title('平面')
r0=abs(z1-z2)<=1;
zz=r0.*z2;yy=r0.*y;xx=r0.*x;
subplot(1,3,3),plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'x')
title('交线')
例4-34利用rotate函数,从不同角度查看函数
t=-2:
.1:
2;
[x,y]=meshgrid(t);
z=x.*exp(-x.^2-y.^2);
subplot(121),mesh(x,y,z);
subplot(122),h=mesh(x,y,z)%返回图形对象的句柄
rotate(h,[-2,-2,0],30,[2,2,0])
subplot(121),surf(ones(10,10))
%subplot(122),h=surf(ones(10,10));rotate(h,[001],45,[100]),rotate3d
fori=1:
45
subplot(122),h=surf(ones(10,10));rotate(h,[001],i,[100]),rotate3d
end
线性代数,数值积分
例1:
求方程组的解
解:
A=[2,3;1,-1];
b=[4;1]
X=inv(A)*b
ans
X=
1.4000
0.4000
方程的解是:
x=1.4,y=0.4
例5-21求齐次线性方程组的通解
解:
Matlab命令为
A=[1-8102;245-1;386-2];系数矩阵
rref(A)行的最简形式
ans=
1040
01-3/4-1/4
0000
将0=0的一行去掉,则原方程组等价于
取
得
取
得
基础解系为
所以方程的通解为
其中k1,k2是任意实数
例5-22求非齐次线性方程组的通解
解:
MATLAB命令为:
B=[1-1-110;1-11-31;1-1-23-1/2];
rref(B)
ans=
1-10-11/2
001-21/2
00000
分析:
原方程组对应的同解方程组为:
其中一个特解为:
再求解对应的齐次线性方程组
可得到一个基础解系:
所以方程的通解为
例求非齐次线性方程组的通解
解:
Matlab命令为
B=[42-12;3-1210;11308];
rref(B)
ans=
结果分析:
行最简形式中最后一行出现了零等于
非零的情况,故方程组无解。
例1计算定积分
与精确值pi相比较
解MATLAB命令为:
h=0.01;x=0:
h:
1;
y=4./(1+x.^2);
formatlong
z1=sum(y(1:
length(x)-1))*h
z2=sum(y(2:
length(x)))*h
formatshort
u1=z1-pi,u2=z2-pi
输出结果:
z1=3.151********3129
z2=3.131********3129
u1=0.0100
u2=-0.0100
例3用三种方法计算定积分
解:
编程如下:
x=0:
0.01:
1;
y=sin(x.^2)./(x+1);
s1=sum(y(1:
100))*0.01
s2=sum(y(2:
101))*0.01
s3=trapz(x,y)
s4=quad('sin(x.^2)./(x+1)',0,1)
运行结果为:
s1=0.1787
s2=0.1829
s3=0.1808
s4=0.1808
按左矩形公式计算结果是0.1787,按右矩形公式计算结果是0.1829,按梯形法和辛普生法计算结果都是0.1808.
建立函数文件jifen.m:
functiony=jifen(x)
y=sin(x.^2)./(x+1);
编程如下:
s=quad('jifen',0,1)
拟合与差值
例6-19:
用以上4种方法对y=cosx在[0,6]上的一维插值效果进行比较。
x=0:
6;
y=cos(x);
xi=0:
.25:
6;
yi1=interp1(x,y,xi,'*nearest');
yi2=interp1(x,y,xi,'*linear');
yi3=interp1(x,y,xi,'*spline');
yi4=interp1(x,y,xi,'*cubic');
plot(x,y,'ro',xi,yi1,'--',xi,yi2,'-',xi,yi3,'k.-',xi,yi4,'m:
')
legend(‘原始数据’,‘最近点插值’,‘线性插值’,’样条插值‘,’立方插值’)
例6-21:
用以上4种方法对
在[-2,2]上的二维多项式插值效果进行比较。
[x,y]=meshgrid(-2:
.5:
2);
z=x.*exp(-x.^2-y.^2);
[x1,y1]=meshgrid(-2:
.1:
2);
z1=x1.*exp(-x1.^2-y1.^2);
figure
(1)
subplot(1,2,1),mesh(x,y,z),title(‘数据点')
subplot(1,2,2),mesh(x1,y1,z1),title(‘函数图象')
[xi,yi]=meshgrid(-2:
.125:
2);
zi1=interp2(x,y,z,xi,yi,'*nearest');
zi2=interp2(x,y,z,xi,yi,'*linear');
zi3=interp2(x,y,z,xi,yi,'*spline');
zi4=interp2(x,y,z,xi,yi,'*cubic');
figure
(2)
subplot(221),mesh(xi,yi,zi1),title(‘最近点插值')
subplot(222),mesh(xi,yi,zi2),title(‘线性插值')
subplot(223),mesh(xi,yi,zi3),title(‘样条插值')
subplot(224),mesh(xi,yi,zi4),title(‘立方插值')
常微分和概率
例7-45:
用数值积分的方法求解微分方程
设初始时间
终止时间
初始条件
分析:
求解
令:
(化为一阶微分方程)即原微分方程化为:
写成矩阵形式为
(化为一阶微分方程)
编写函数文件exf.m
functionxdot=exf(t,x)
u=1-(t.^2)/(2*pi);
xdot=[0,1;-1,0]*x+[01]'*u;
定义另外一个函数为主函数
clf;
t0=0;
tf=3*pi;
x0t=[0;0];
[t,x]=ode23('exf',[t0,tf],x0t)
y=x(:
1);%[t,x]中求出的x是按列排列,故用ode23求出x后只要第
一列即为y
y2=-1/2*(-2*pi-2+t.^2)/pi-(pi+1)/pi*cos(t);
clf,
plot(t,y,'-',t,y2,'o')
在100人的团体中,如果不考虑年龄的差异,研究是否有两个以上的人生日相同。
假设每人的生日在一年365天中的任意一天是等可能的,那么随机找个人(<365)
问这些人生日各不相同的概率是多少?
至少有两个人生日相同的概率为多少?
forn=1:
100
p0(n)=prod(365:
-1:
365-n+1)/365^n;
p1(n)=1-p0(n);
end
n=1:
100;
plot(n,p0,n,p1,'--')
xlabel(‘人数’),ylabel(‘概率’)
legend(‘生日各不相同的概率’,‘至少两人相同的概率’)
axis([0100-0.11.1]),gridon
例8-19某班(共有120名学生)的高等数学成绩如下:
746378768956709789947688
658372413972736814764570
904654617576495778666474
788786734767216679676865
568466736872766570945365
777853745950986789786392
548784806364856669696054
753330627465847355857576
817183725684767567653594
594745677536788294708475
根据以上数据作出该门课程成绩的频数表和直方图。
解:
(1)数据输入:
v方法1:
在Matlab的交互环境下直接输入;
v方法2:
将以上数据以一列的形式存为A.txt文件,用
(2)用hist命令作频数表和直方图:
(区间个数为5,可省略)
[N,X]=hist(A,5)120名学生高数成绩的频数表;
hist(A,5)120名学生高数成绩的直方图;
matleb命令:
loadA.txt↙
disp(‘高数成绩的频数表’),[N,X]=hist(A,5)%N为频数
hist(A,5)%直方图
[N,X]=hist(A,5)
N=
310226025
X=
22.400039.200056.000072.800089.6000
a1=min(A);a2=max(A);
disp([‘成绩最小值’,blanks(4),‘最大值'])
disp([a1,a2])
成绩最小值最大值
1498
例8-20求例8-19中A的均值、中位数、极差、方差
和标准差。
解:
在命令窗口输入:
M=[mean(A),median(A),range(A),var(A),std(A)]
M=68.958371.500084.0000249.569715.7978
课后习题
p16习题4
1.程序如下:
functionfun(n)
clc
sum=0;
s=1;
fori=1:
n
s=s*i;
sum=sum+s;
end
sum
结果如下:
sum=
2.5613e+018
p27习题3
程序如下:
A=[52;91]
A=
52
91
>>B=[12;92]
B=
12
92
结果:
>>A>B
ans=
10
00
>>A==B
ans=
01
10
>>A
ans=
00
01
>>(A==B)&(A
ans=
00
00
>>(A==B)&(A>B)
ans=
00
00
p34习题1
程序如下:
>>s=0;
i=1;
whileabs(1/(2*i-1))>=exp(-6)
s=s+1/(2*i-1)*(-1)^(i+1);
i=i+1;
end
p=4*s;
结果如下:
>>p
p=
3.1366
p79习题1
程序如下:
clc
clf
x=0:
pi/50:
4*pi;
y1=exp(x/3).*sin(3*x);
y2=exp(x/3);
y3=-exp(x/3);
plot(x,y1,'b*',x,y2,'r-.',x,y3,'r-.'),gridon
结果如下:
p79习题3
2.程序如下:
clc
x1=-pi:
pi/50:
4*pi;
x2=pi:
pi/40:
4*pi;
x3=1:
0.1:
8;
y1=x1.*cos(x);
y2=x2.*tan(1./x2).*sin(x2.^3);
y3=exp(1./x3).*sin(x3);
subplot(1,3,1),plot(x1,y1)
axis([-pipi-44]),gridon
title('y=xcosx'),gtext('y=xcosx')
legend('y=xcosx')
xlabel('x'),ylabel('y=xcosx'),
subplot(1,3,2),plot(x2,y2)
axis([pi4*pi-44]),gridon
title('y=xtan(1/x)sin(x^3)'),gtext('y=xtan(1/x)sin(x^3)')
legend('y=xtan(1/x)sin(x^3)')
xlabel('x'),ylabel('y=xtan(1/x)sin(x^3)'),
subplot(1,3,3),plot(x3,y3)
axis([18-44]),gridon
title('exp(1/x)sinx'),gtext('exp(1/x)sinx')
legend('exp(1/x)sinx')
xlabel('x'),ylabel('exp(1/x)sinx'),
结果:
p79习题5
3.程序如下:
t=0:
pi/10:
20*pi;
x=t.*cos(pi/6*t);
y=t.*sin(pi/6*t);
z=2*t;
plot3(x,y,z),gridon,title('圆锥螺线')
xlabel('x=tcos(pi*t/6)')
ylabel('y=tsin(pi*t/6)')
zlabel('z=2t')
gtext('圆锥螺线')
结果如下:
p79习题9
程序如下:
t=-2:
0.1:
2;
[x,y]=meshgrid(t);
z1=5-x.^2-y.^2;
subplot(1,3,1)
mesh(x,y,z1)
z2=3*ones(size(x));
subplot(1,3,2),mesh(x,y,z2)
r0=abs(z1-z2)<=0.01;
zz=r0.*z2;yy=r0.*y;xx=r0.*x;
subplot(1,3,3),plot3(xx(r0==1),yy(r0==1),zz(r0==1))
结果:
P114习题12
程序及结果如下:
>>A=[123456;
789101112;
131415161718;
192021222324;
252627282930;
313233343536;]
A=
123456
789101112
131415161718
192021222324
252627282930
313233343536
>>A'
ans=
1713192531
2814202632
3915212733
41016222834
51117232935
61218243036
>>det(A)
ans=
3.5004e-059
>>rank(A)
ans=
2
>>rref(A)
ans=
10-1-2-3-4
012345
000000
000000
000000
100000
习题14
程序如下:
>>A=[211;121;112];
>>p=poly(A)
p=
1.0000-6.00009.0000-4.0000
>>d=eig(A)
d=
1.0000
1.0000
4.0000
>>[V,D]=eig(A)
V=
0.40820.70710.5774
0.4082-0.70710.5774
-0.816500.5774
D=
1.000000
01.00000
104.0000
习题21(1,
解:
matlab命令为:
>>A=[112-1;-1130;2-34-1];
>>rref(A)
ans=
1.000000-0.5600
01.00000-0.2000
001.0000-0.1200
结果分析:
可以看出系数矩阵A的秩为3,小于未知量的个数4,所以有无穷多个解,原方程组对应的同解方程组为:
x1=0.56x4
x2=0.2x4
x3=0.12x4
取x4=1,解的方程组的基础解系为:
0.56
0.2
ξ=0.12
1
所以方程组的通解为:
X10.56
X20.2
X3=k0.12其中k为任意实数
X41
习题21(2
解:
matlab命令为:
>>B=[1-1-110;1-11-31;1-1-23-1/2];
>>rref(B)
ans=
1.0000-1.00000-1.00000.5000
001.0000-2.00000.5000
00000
结果分析:
可以看出增广矩阵的秩为2,等于系数矩阵的秩,而小于未知量的个数4,所以方程组有无穷多个解,原方程组对应的同解方程为:
X1=x2+x4+0.5
X3=2x4+0.5
可以找到其中一个特解为:
2.5
η=1
2.5
1
再求解对应的齐次线性方程组
X1=x2+x4
X3=2x4
可以得到一个基础解系:
11
10
ξ1=0,ξ2=2
01
因此方程组的通解为:
X1112.5
X2101
X3=c10+c22+2.5
X4011其中c1,c2为任意实数
P167习题17
(2)
解:
建立函数文件jifen.m:
functiony=jifen(x)
y=x.*sin(x)./(1+cos(x).*cos(x));
编程如下:
x=0:
pi/100:
pi;
y=x.*sin(x)./(1+cos(x).*cos(x));
t=length(x);
s1=sum(y(1:
(t-1)))*pi/100
s2=sum(y(2:
t))*pi/100
s3=trapz(x,y)
s4=quad('jifen',0,pi)
运行结果为:
s1=
2.4673
s2=
2.4673
s3=
2.4673
s4=
2.4674
习题18
解:
matlab命令为:
h=pi/1000;
x=0:
h:
pi/4;
y=1./(1-sin(x));
formatlong
t=length(x);
z1=sum(y(1:
(t-1)))*h
z2=sum(y(2:
t))*h
z3=trapz(x,y)
z4=quad('1./(1-sin(x))',0,pi/4)
formatshort
u1=z1-sqrt
(2),u2=z2-sqrt
(2),u3=z3-sqrt
(2),u4=z4-sqrt
(2)
输出结果为:
z1=
1.410427281389369
z2=
1.418011756981117
z3=
1.4142195191