北邮计算物理实践ch3数组作业

北邮计算物理实践ch3数组作业

1.

d=pi/100;
t=0:d:2*pi;
x=sin(t);
dt=10*eps;                      %
x_eps=sin(t+dt);
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')  

t1n = ndims(t1)
t1s = size(t1)
t1l = length(t1)  

Tn = ndims(T)
Ts = size(T)
Tl = length(T)  

Fn = ndims(F)
Fs = size(F)
Fl =length(F)  

T12 = (t1==t2)  

T12 =
  1×10 logical 数组
   1   1   1   1   1   1   0   1   1   1  
T12(7)=0,说明t1在第七个元素处不等于t2，我们可以来看一下：
format long  
t1_7 = t1(7)  

t1_7 =
   2.094395102393196  
t2_7 = t2(7)  

t2_7 =
   2.094395102393195  
可以看出t1与t2的第七个元素在最后一位不同

2.

A = reshape(1:18,3,6)  

A =
     1     4     7    10    13    16
     2     5     8    11    14    17
     3     6     9    12    15    18  
A([2,4,8,16]) = NaN  

A =
     1   NaN     7    10    13   NaN
   NaN     5   NaN    11    14    17
     3     6     9    12    15    18  
A(:,[4,5]) = Inf  

A =
     1   NaN     7   Inf   Inf   NaN
   NaN     5   NaN   Inf   Inf    17
     3     6     9   Inf   Inf    18  

3.

rng('default')
A = rand(3,5)  

A =
    0.8147    0.9134    0.2785    0.9649    0.9572
    0.9058    0.6324    0.5469    0.1576    0.4854
    0.1270    0.0975    0.9575    0.9706    0.8003  

[ar,aj] = find(A>0.5); 
display(['全下标索引的行号',int2str(ar')])
display(['全下标索引的列号',int2str(aj')])  

全下标索引的行号1  2  1  2  2  3  1  3  1  3
全下标索引的列号1  1  2  2  3  3  4  4  5  5  

ak = find(A>0.5);
display(['单下标',int2str(ak')])  

单下标1   2   4   5   8   9  10  12  13  15  

4.

fprintf('标量循环法的纵轴是y1,数组运算法是y2');
A = linspace(0,10,100);
N = length(A);
B = zeros(size(A));
for k = 1:N
    B(k) = 1-exp(-0.5*A(k))*cos(2*A(k));
end
subplot(1,2,1)
plot(A,B)
xlabel('t');
ylabel('y1');
grid on

C  = 1-exp(-0.5*A).*cos(2*A);
subplot(1,2,2)
plot(A,C)
xlabel('t');
ylabel('y2');
grid on  

标量循环法的纵轴是y1,数组运算法是y2

5.

A = magic(3); B = rand(3);  
(1)
AB1 = A.*B == B.*A  

AB1 =
  3×3 logical 数组
   1   1   1
   1   1   1
   1   1   1  
相同

(2)
AB2 = A*B == A.*B  

AB2 =
  3×3 logical 数组
   0   0   0
   0   0   0
   0   0   0  

不相同，因为A*B是正常代数上的矩阵乘法，而点乘是对应元素相乘

(3)
AB3 = A*B == B*A  

AB3 =
  3×3 logical 数组
   0   0   0
   0   0   0
   0   0   0  
不同，显然代数中的A*B矩阵和B*A矩阵不同

(4)
AB4 = A.\B == B./A  

AB4 =
  3×3 logical 数组
   1   1   1
   1   1   1
   1   1   1  
相同，因为点除都是B矩阵元素除以A矩阵对应元素

(5)
AB5 = A\B == B/A  

AB5 =
  3×3 logical 数组
   0   0   0
   0   0   0
   0   0   0  
不同，因为A\B是Ax = B中的解x，而B/A是xA = B中的解，在第(3)题中已经论证Ax与xA结果不同

(6)
AB6 = A*A\B-B == A*(A\B)-B  

AB6 =
  3×3 logical 数组
   0   0   0
   0   0   0
   0   0   0  
不同，第二个接近于0

(7)
AB7 = A\eye(3) == eye(3)/A  

AB7 =
  3×3 logical 数组
   1   1   0
   0   1   1
   0   1   0  

运算结果相同，因为得到的都是A的逆矩阵

6.

A = [1 2; 3 4]  

A =
     1     2
     3     4  
(1)
B1 = A.^0.5  

B1 =
    1.0000    1.4142
1.7321    2.0000  

B2 = 0.5.^A  

B2 =
    0.5000    0.2500
    0.1250    0.0625  

B3 = A^0.5  

B3 =
   0.5537 + 0.4644i   0.8070 - 0.2124i
   1.2104 - 0.3186i   1.7641 + 0.1458i  

B4 = 0.5^A  

B4 =
    0.9910   -0.4422
   -0.6634    0.3276  

(2)
A1 = B1.^2;  
A2 = log(B2)/log(0.5);  
A3 = B3^2;  
A4 = logm(B4)/logm(0.5);  

(3)
其实这些矩阵实部都一样，虚部为0而已

7.

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

(1)
Z中有91+90 = 181 个非数数据，
L = find(isnan(Z));     %NaN ~= NaN,所以不能用find(Z == NaN)
n = length(L)  

n =
   181 

surf(X,Y,Z)
shading interp  

(2)
x = -3*pi:pi/15:3*pi;  
y = x;  
[X,Y] = meshgrid(x,y); 
X2 = X + find(X == 0)*realmin;
Y2 = Y + find(Y == 0)*realmin; 
warning off;  
Z2 = sin(X2).*sin(Y2)./X2./Y2; 
surf(X2,Y2,Z2) 
shading interp  

8.

(1)标量循环+条件分支法：
x = -1.5:0.001:1.5; y = -3:0.001:3;  
[X,Y] = meshgrid(x,y);  
xl = length(x); yl = length(y);
z = zeros(xl,yl);
for ii = 1:xl
    for jj = 1:yl
        if (    x(ii)+y(jj)<=-1)
            z(ii,jj) = 0.546*exp(-0.75*y(jj)^2-3.75*x(ii)^2+1.5*x(ii));
        elseif (x(ii)+y(jj)>-1 & x(ii)+y(jj)<=1)
            z(ii,jj) = 0.758*exp(-y(jj)^2-6*x(ii)^2);
        else
            z(ii,jj) = 0.546*exp(-0.75*y(jj)^2-3.75*x(ii)^2-1.5*x(ii));
        end
    end
end

surf(X,Y,z')
shading interp  

(2)数组混合
x = -1.5:0.001:1.5; y = -3:0.001:3;
[X,Y] = meshgrid(x,y);

Z1 = 0.546*exp(-0.75*Y.^2-3.75*X.^2+1.5*X).*(X+Y<=-1);
Z2 = 0.758*exp(-Y.^2-6*X.^2).*(X+Y>-1 & X+Y<=1);
Z3 = 0.546*exp(-0.75*Y.^2-3.75*X.^2-1.5*X).*(X+Y>1);

Z = Z1+Z2+Z3; 

surf(X,Y,Z)
shading interp  

(3)比较两种方法的元素相对误差
format long  
m1 = abs(z'-Z);
m2 = max(m1(:))  

9.

(1)
function D = matrix_old_multiply(A,B)
% 矩阵的乘法
%   对应行列相乘相加
[a,b] = size(A);
[c,d] = size(B);
if b~=c
    error('A矩阵的列数不等于B矩阵的行数，不能相乘啦！')
end

for ii = 1:a
    for jj = 1:d
        new_i = 0; new_j = 0;
        for k = 1:b
            new_i = new_i + real(A(ii,k))*real(B(k,jj))...
                -imag(A(ii,k))*imag(B(k,jj));
            new_j = new_j + real(A(ii,k))*imag(B(k,jj))...
                +imag(A(ii,k))*real(B(k,jj));
        end
        D(ii,jj) = new_i + i*new_j;
    end
end

end

计算G = (AB-C)*D只需调用两次上面的函数即可
rng default
A=randn(50,70)+1i*randn(50,70);
B=randn(70,60)+1i*randn(70,60);
C=randn(50,60)+1i*randn(50,60);
D=randn(60,1)+1i*randn(60,1);  

K = matrix_old_multiply(A,B);
G = matrix_old_multiply(K-C,D);  

(2)

G2 = (A*B-C)*D;
N1 = norm(real(G2))
N2 = norm(imag(G2))  

N1 =
     6.352169123696923e+02
N2 =
     8.048619385149246e+02