A_sinA =
0.8415 0.9093 0.1411 -0.7568 -0.9589 -0.2794 0.6570 0.9894 0.4121 A_sinM =
-0.6928 -0.2306 0.2316 -0.1724 -0.1434 -0.1143
0.3479 -0.0561 -0.4602
4.4 奇异值分解
4.4.1 奇异值分解和矩阵结构 4.4.1.1 奇异值分解定义
4.4.1.2 矩阵结构的奇异值分解描述 4.4.2 线性二乘问题的解 4.4.2.1 矩阵除运算的广义化 4.4.2.2 线性模型的最小二乘解
【例4.4.2.2-1】对于超定方程y?Ax,进行三种解法比较。其中A取MATLAB库中的特殊函数生成。 (1)
A=gallery(5);A(:,1)=[];y=[1.7 7.5 6.3 0.83 -0.082]';
x=inv(A'*A)*A'*y,xx=pinv(A)*y,xxx=A\\y
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 5.405078e-018. x =
3.4811e+000 5.1595e+000 9.5340e-001 -4.6569e-002 xx =
3.4759e+000 5.1948e+000 7.1207e-001 -1.1007e-001
Warning: Rank deficient, rank = 3 tol = 1.0829e-010. xxx =
3.4605e+000 5.2987e+000 0 -2.9742e-001
(2)
nx=norm(x),nxx=norm(xx),nxxx=norm(xxx) nx =
6.2968e+000 nxx =
6.2918e+000
6
nxxx =
6.3356e+000
(3)
e=norm(y-A*x),ee=norm(y-A*xx),eee=norm(y-A*xxx) e =
6.9863e-001 ee =
4.7424e-002 eee =
4.7424e-002
4.5 函数的数值导数和切平面
4.5.1 法线
【例4.5.1-1】曲面法线演示。
y=-1:0.1:1;x=2*cos(asin(y)); [X,Y,Z]=cylinder(x,20);
surfnorm(X(:,11:21),Y(:,11:21),Z(:,11:21)); view([120,18]) 1.210.80.60.40.20-0.2-3-5-20-1015图 4.5.1-1 4.5.2 偏导数和梯度 4.5.2.1 理论定义 4.5.2.2 数值计算指令
【例4.5.2.2-1】用一个简单矩阵表现diff和gradient指令计算方式。
F=[1,2,3;4,5,6;7,8,9] Dx=diff(F)
Dx_2=diff(F,1,2) [FX,FY]=gradient(F)
[FX_2,FY_2]=gradient(F,0.5) F =
1 2 3 4 5 6 7 8 9 Dx =
3 3 3 3 3 3 Dx_2 =
1 1
7
1 1 1 1 FX =
1 1 1 1 1 1 1 1 1 FY =
3 3 3 3 3 3 3 3 3 FX_2 =
2 2 2 2 2 2 2 2 2 FY_2 =
6 6 6 6 6 6 6 6 6
【例 4.5.2.2-2】研究偶极子(Dipole)的电势(Electric potential)和电场强度(Electric field density)。设在(a,b)处有电荷?q,在(?a,?b)处有电荷?q。那么在电荷所在平面上任
?11何一点的电势和场强分别为V(x,y)?(?),E???V。其中
4??0r?r?1r??(x?a)2?(y?b)2,r??(x?a)2?(y?b)2。?9?109。又设电荷
4??0q?2?10?6,a?1.5,b??1.5。
qclear;clf;q=2e-6;k=9e9;a=1.5;b=-1.5;x=-6:0.6:6;y=x; [X,Y]=meshgrid(x,y);
rp=sqrt((X-a).^2+(Y-b).^2);rm=sqrt((X+a).^2+(Y+b).^2); V=q*k*(1./rp-1./rm); [Ex,Ey]=gradient(-V);
AE=sqrt(Ex.^2+Ey.^2);Ex=Ex./AE;Ey=Ey./AE; cv=linspace(min(min(V)),max(max(V)),49); contourf(X,Y,V,cv,'k-') %axis('square')
title('\\fontname{隶书}\\fontsize{22}偶极子的场'),hold on quiver(X,Y,Ex,Ey,0.7) plot(a,b,'wo',a,b,'w+')
plot(-a,-b,'wo',-a,-b,'w-')
xlabel('x');ylabel('y'),hold off ???×??6420y-2-4-6-6-4-20x246图 4.5.2.2-1 4.6 函数的零点
8
4.6.1 多项式的根 4.6.2 一元函数的零点
4.6.2.1 利用MATLAB作图指令获取初步近似解 4.6.2.2 任意一元函数零点的精确解
【例 4.6.2.2-1】通过求(1)
y=inline('sin(t)^2*exp(-a*t)-b*abs(t)','t','a','b');
%<1>
f(t)?(sin2t)e?at?bt的零点,综合叙述相关指令的用法。
(2)
a=0.1;b=0.5;t=-10:0.01:10;
y_char=vectorize(y); % <3> Y=feval(y_char,t,a,b);
clf,plot(t,Y,'r');hold on,plot(t,zeros(size(t)),'k'); xlabel('t');ylabel('y(t)'),hold off 10-1y(t)-2-3-4-5-10-50t510图4.6-1
(3) 由于Notebook中无法实现zoom、ginput指令涉及的图形和鼠标交互操作,因此下面指令必须在MATLAB指令窗中运行,并得到如图4.6-2所示的局部放大图及鼠标操作线。
zoom on
[tt,yy]=ginput(5);zoom off
图 4.6-2
tt tt =
-2.0032 -0.5415 -0.0072
9
0.5876 1.6561
(4)
[t4,y4,exitflag]=fzero(y,tt(4),[],a,b) t4 =
0.5993 y4 = 0
exitflag = 1
%<11>
(5)
[t3,y3,exitflag]=fzero(y,tt(3),[],a,b) t3 =
-0.5198 y3 =
-5.5511e-017 exitflag = 1
(6)
op=optimset('fzero') op =
ActiveConstrTol: [] ......
Display: 'notify' ......
TolX: 2.2204e-016 TypicalX: []
op=optimset('tolx',0.01); op.TolX ans = 0.0100
(7)
[t4n,y4n,exitflag]=fzero(y,tt(4),op,a,b) t4n = 0.6042 y4n = 0.0017 exitflag = 1
4.6.3 多元函数的零点
【例 4.6.3-1】求解二元函数方程组?(1)
x=-2:0.5:2;y=x;[X,Y]=meshgrid(x,y); F1=sin(X-Y);F2=cos(X+Y); v=[-0.2, 0, 0.2]; contour(X,Y,F1,v)
hold on,contour(X,Y,F2,v),hold off
?f1(x,y)?sin(x?y)?0的零点。
?f2(x,y)?cos(x?y)?0 10