matlab数值计算 1

第三章课程复习。

1．要求在闭区间上产生具有10个等距采样点的一维数组。试用两种不同的指令实现。

linspace(0,2*pi,10)

0:2*pi/9:2*pi

2．在使用123作为rand随机数发生器的初始化状态的情况下，写出产生长度为1000的“等概率双位（即取-1，+1）取值的随机码”程序指令，并给出 -1码的数目。

答案〗a=randsrc(1000,1,[-1,1],123);

b=sum(a==-1)

3．已知矩阵，运行指令b1=a.^(0.5), b2=a^(0.

5), 可以观察到不同运算方法所得结果不同。（1）请分别写出根据b1, b2恢复原矩阵a的程序。（2）用指令检验所得的两个恢复矩阵是否相等。

a=[1 2;3 4]

b1=a.^(0.5), b2=a^(0.5)

b1.^2b2^2

4．先运行clear,format long,rand('state',1),a=rand(3,3)，然后根据a写出两个矩阵：一个对角阵b，其相应元素由a的对角元素构成；另一个矩阵c，其对角元素全为0，而其余元素与对应的a阵元素相同。

clear,format long,rand('state',1),a=rand(3,3)

b=diag(diag(a))

c=a-b

5．先运行指令x=-3*pi:pi/15:3*pi; y=x; [x,y]=meshgrid(x,y); warning off; z=sin(x).

*sin(y)./x./y; 产生矩阵z。

（1）请问矩阵z中有多少个“非数”数据？（2）用指令surf(x,y,z); shading interp观察所绘的图形。（3）请写出绘制相应的“无裂缝”图形的全部指令。

x=-3*pi:pi/15:3*pi; y=x; [x,y]=meshgrid(x,y); warning off; z=sin(x).*sin(y)./x./y;

a=sum(sum(isnan(z)))

surf(x,y,z); shading interp

xx=x+(x==0)*eps;

yy=y+(y==0)*eps;

x,y]=meshgrid(xx,yy);

z=sin(x).*sin(y)./x./y;

surf(x,y,z); shading interp;

答案〗6．下面有一段程序，企图用来解决如下计算任务：有矩阵，当依次取10, 9, 8, 7, 6, 5, 4, 3, 2, 1时，计算矩阵“各列元素的和”，并把此求和结果存放为矩阵sa的第k行。例如时，a阵为，此时它各列元素的和是一个行数组，并把它保存为sa的第3行。

问题：该段程序的计算结果对吗？假如计算结果不正确，请指出错误发生的根源，并改正之。

for k=10:-1:1

a=reshape(1:10*k,k,10);

sa(k,:)sum(a);

endsa

与符号计算相比，数值计算在科研和工程中的应用更为广泛。

和十年、二十年前相比较，在当今计算机软硬件的支持下，人们所能拥有的计算能力已经得到了巨大的提升。。

matlab数值计算中，没有专门的求极限和导数的指令。

原因：数值精度有限，不能表示无穷小量，不能准确描述一个数的邻域。

例4.1-1】设，，试用机器零阈值eps替**论0计算极限，。

x=eps;

l1=(1-cos(2*x))/x*sin(x)),

l2=sin(x)/x,

syms t

f1=(1-cos(2*t))/t*sin(t));

f2=sin(t)/t;

ls1=limit(f1,t,0)

ls2=limit(f2,t,0)

cleark=logspace(0,14,15);

x=k*eps;

f1=(1-cos(2.*x)).x./sin(x)/2;

f2=sin(x)./x;

format short e

disp([blanks(7),'x',blanks(11),'f1',blanks(11),'f2'])

disp([x',f1',f2

clfn1=sum(f1<10^-10

subplot(2,1,1)

semilogx(x(1:n1),f1(1:n1),'ro','linewidth',3

hold on

semilogx(x(n1+1:end),f1(n1+1:end),'b.-'linewidth',3)

hold off

text(10^-14,0.2,'错误区段'),text(10^-6,0.8,'正确区段')

axis([10^-20,1,0,1.1])

grid on

title('f2')

subplot(2,1,2)

semilogx(x,f2,'b.-'linewidth',3)

axis([10^-20,1,0,1.1])

title('f1')

xlabel('x')

grid on

shg 理论分析表明：，

借助符号计算所求的极限与理论值一致；

用数值法近似计算的极限与理论不一致。

提醒：除非数值近似法求的极限经过理论验证，否则绝不要借助数值法求取极限。

函数在点处的导数。

例4.1-2】已知，求该函数在区间中的近似导函数。

d=pi/100;

t=0:d:2*pi;

x=sin(t);

dt=5*eps

x_eps=sin(t+dt自变量的增量为5eps。

dxdt_eps=(x_eps-x)/dt近似数值导数。

plot(t,x,'linewidth',5)

hold on

plot(t,dxdt_eps)

hold off

legend('x(t)',dx/dt')

xlabel('t')

4.1-1 增量过小引起有效数字严重丢失后的毛刺曲线。

x_d=sin(t+d

dxdt_d=(x_d-x)/d以d=pi/100为增量算得的数值导数。

plot(t,x,'linewidth',5)

hold on

plot(t,dxdt_d)

hold off

legend('x(t)',dx/dt')

xlabel('t')

图 4.1-2 增量适当所得导函数比较光滑。

结论： matlab提供的与“求导”概念有关的指令：

dx = diff(x求差分。

当x是向量时，dx = x(2:n) –x(1:n-1)；

当x是矩阵时，dx = x(2:n,:)x(1:n-1,:)

dx的长度比x的“长度”短一个元素（或一行）。

x = 1 2 3 4 5];

y = diff(x)

fx = gradient(f求梯度。

fx,fy] =gradient(f) 二元函数的梯度。

当f是向量时，fx(1) =f(2) –f(1),

fx(end) =f(end) –f(end-1),

fx(2:end-1) =f(3:end) –f(1:end-2))/2,

fx的“长度”与f相同。

当f是矩阵时，fx和fy是与f同样大小的矩阵。

fx的每行给出f相应行元素间的“梯度”；

fy的每列给出f相应列元素间的“梯度”。

例4.1-3】已知，采用diff和gradient计算该函数在区间中的近似导函数。

clfd=pi/100

t=0:d:2*pi;

x=sin(t);

dxdt_diff=diff(x)/d;

dxdt_grad=gradient(x)/d;

subplot(1,2,1)

plot(t,x,'b')

hold on

plot(t,dxdt_grad,'m','linewidth',8)

plot(t(1:end-1),dxdt_diff,'.k','markersize',8)

axis([0,2*pi,-1.1,1.1])

title('[0, 2\pi]')

legend('x(t)',dxdt_',dxdt_',location','north')

xlabel('t'),box off

hold off

subplot(1,2,2)

kk=(length(t)-10):length(t); 最后11个数据的下标。

hold on

plot(t(kk),dxdt_grad(kk),'om','markersize',8) %微观上看，不仅数值上有差异，plot(t(kk-1),dxdt_diff(kk-1),'k','markersize',8) %而且diff没有给出最后一点的导数。

title('[end-10, end]')

legend('dxdt_',dxdt_',location','southeast')

xlabel('t'),box off

hold off

图 4.1-3 diff和gradient求数值近似导数的异同比较。

sx=sum(x沿列方向求和。

scs=cumsum(x沿列方向求累计和。

st=trapz(x,y采用梯形法沿列方向求函数y关于自变量x的积分。

sct=cumtrapz(x,y) 采用梯形法沿列方向求函数y关于自变量x的累计积分。

matlab数值计算 1

3MATLAB数值计算

1MATLAB数值运算

数值分析matlab版第二章

其他用户还读了