Python和Matlab之间的矩阵计算不同答案

【问题标题】：Matrix Calculation Different between Python and MatlabPython和Matlab之间的矩阵计算不同
【发布时间】：2020-07-24 08:10:10
【问题描述】：

我正在将有限元分析代码从 Matlab 转移到 Python。

当我尝试求解位移 U = F/K 时，在最后一步出现问题。我检查了 Matlab 和 Python 计算的 F 和 K 是否相同。但是得到的U是不同的。当矩阵大小达到19684*19684时，重要的数字会完全不同，甚至符号会相反。

我的问题是针对相同的 K 和 F 矩阵，为什么 Python 和 MATLAB 会导致 U = F/K 的答案不同。请注意，如果 K 和 F 的大小很小（例如 1000* 1000,1000* 1），则所得 U = F/K 与使用 Python 和 MATLAB 非常接近。但是，如果这两个矩阵的大小更大，那么得到的 U 就会完全不同。

谁能建议这里有什么问题，我该如何解决？

下面是 Matlab 和 Python 的代码。

Matlab

clear all; clc; close all
nelx=80; 
nely=80; 
beta = 1;
eta = 0.47;
volrate=0.20; 
E0 = 1;
Emin = 0.000001;
nu = 0.3;
a = 0.005;
b = 0.005;
h = 0.01;
q = -0.001;
penal = 3;
[KE1]= ElementStiffness(E0,nu,a,b,h);
TM3 = zeros(nely,nelx);
TM4 = zeros(nely+1,nelx+1);
TM3 = TM3+0.01;
% TM4([98:101 199:202 300:303 401:404]) = 0.01;
TM4([1]) = 0.01;
T2 = zeros((nely+1)*(nelx+1),1); 
num = 0;
for i=1:nely+1
    for j =1:nelx+1
        num = num+1;
        T2(num,1) = TM4(j,i);
    end
end
xTilde = repmat(volrate,[nely nelx]);
xPhys = ((tanh(beta*eta)+tanh(beta*(xTilde-eta)))/(tanh(beta*eta)+tanh(beta*(1-eta))));


K=sparse(3*(nelx+1)*(nely+1),3*(nelx+1)*(nely+1)); 
F=sparse(3*(nely+1)*(nelx+1),1);
U1=sparse(3*(nely+1)*(nelx+1),1); 

Fe=q*a*b*[1 b/3 -a/3 1 b/3 a/3 1 -b/3 a/3 1 -b/3 -a/3]'; 

for ely=1:nely 
    for elx=1:nelx 
        n1 = (nely+1)*(elx-1)+ely;
        n2 = (nely+1)* elx +ely;
        edof = [3*n1+1; 3*n1+2; 3*n1+3;3*n2+1;3*n2+2;3*n2+3;3*n2-2;3*n2-1; 3*n2; 3*n1-2;3*n1-1; 3*n1];
        K(edof,edof)=K(edof,edof)+(Emin+(1-Emin)*xPhys(ely,elx)^penal)*KE1; 
        if TM3(ely,elx)>0 
            F(edof)=F(edof)+Fe;
        end
    end 
end 


dd = find(T2 == 0.01); %T2=0-fixed node
fixeddofs = [(dd-1)*3+1 (dd-1)*3+2 (dd-1)*3+3]; 
alldofs   =1:3*(nely+1)*(nelx+1); 
freedofs  =setdiff(alldofs,fixeddofs); 

U1(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:); 
U1(fixeddofs,:)=0;
U1matlab = full(U1);

function [KE1]= ElementStiffness(E,nu,a,b,h)
syms x y 
N1  = (1 - x/a)*(1 - y/b)*(2 - x/a - y/b - (x/a)^2 - (y/b)^2)/8;
Nx1 = b*(1 - x/a)*(1 - y/b)*(1 - (y/b)^2)/8;
Ny1 =  - a*(1 - x/a)*(1 - y/b)*(1 - (x/a)^2)/8;
N2  = (1 + x/a)*(1 - y/b)*(2 + x/a - y/b - (x/a)^2 - (y/b)^2)/8;
Nx2 = b*(1 + x/a)*(1 - y/b)*(1 - (y/b)^2)/8;
Ny2 =   a*(1 + x/a)*(1 - y/b)*(1 - (x/a)^2)/8;
N3  = (1 + x/a)*(1 + y/b)*(2 + x/a + y/b - (x/a)^2 - (y/b)^2)/8;
Nx3 = -b*(1 + x/a)*(1 + y/b)*(1 - (y/b)^2)/8;
Ny3 =   a*(1 + x/a)*(1 + y/b)*(1 - (x/a)^2)/8;
N4  = (1 - x/a)*(1 + y/b)*(2 - x/a + y/b - (x/a)^2 - (y/b)^2)/8;
Nx4 = -b*(1 - x/a)*(1 + y/b)*(1 - (y/b)^2)/8;
Ny4 =  - a*(1 - x/a)*(1 + y/b)*(1 - (x/a)^2)/8;
N=[N1 Nx1 Ny1 N2 Nx2 Ny2 N3 Nx3 Ny3 N4 Nx4 Ny4];
B1=diff(N,x,2);
B2=diff(N,y,2);
B3=diff(N,x,y);
B=-[B1;B2;2*B3];
D = (E/(1-nu*nu))*[1, nu, 0 ; nu, 1, 0 ; 0, 0, (1-nu)/2];
BD = transpose(B)*D*B;
KE1=int(int(BD,x,-a,a),y,-b,b)*h^3/12;
KE1=double(KE1);
end

Python

from __future__ import division 
import numpy as np 
from scipy.sparse.linalg import spsolve
from scipy import sparse

nelx=80
nely=80
VF = 0.2
beta = 1
eta = 0.47
penal = 3 
a = 0.005
b = 0.005
q = -0.001
Emin = 1e-6
TM3 = np.zeros((nely,nelx))
TM3 = TM3+0.01
TM4 = np.zeros((nely+1,nelx+1))
TM4[0,0] = 0.01
T2 = np.zeros(((nely+1)*(nelx+1),1))
num = 0
for i in range(nely+1):
    for j in range(nelx+1):       
        T2[num,0] = TM4[j,i]
        num = num+1
KE1 = [[ 9.67032967e-03,  2.23443223e-05, -2.23443223e-05,
    -4.17582418e-03,  5.12820513e-06, -1.95970696e-05,
    -1.31868132e-03,  7.87545788e-06, -7.87545788e-06,
    -4.17582418e-03,  1.95970696e-05, -5.12820513e-06],
   [ 2.23443223e-05,  1.39194139e-07, -2.74725275e-08,
     5.12820513e-06,  4.39560440e-08,  0.00000000e+00,
    -7.87545788e-06,  3.47985348e-08,  0.00000000e+00,
    -1.95970696e-05,  5.67765568e-08,  0.00000000e+00],
   [-2.23443223e-05, -2.74725275e-08,  1.39194139e-07,
     1.95970696e-05,  0.00000000e+00,  5.67765568e-08,
     7.87545788e-06,  0.00000000e+00,  3.47985348e-08,
    -5.12820513e-06,  0.00000000e+00,  4.39560440e-08],
   [-4.17582418e-03,  5.12820513e-06,  1.95970696e-05,
     9.67032967e-03,  2.23443223e-05,  2.23443223e-05,
    -4.17582418e-03,  1.95970696e-05,  5.12820513e-06,
    -1.31868132e-03,  7.87545788e-06,  7.87545788e-06],
   [ 5.12820513e-06,  4.39560440e-08,  0.00000000e+00,
     2.23443223e-05,  1.39194139e-07,  2.74725275e-08,
    -1.95970696e-05,  5.67765568e-08,  0.00000000e+00,
    -7.87545788e-06,  3.47985348e-08,  0.00000000e+00],
   [-1.95970696e-05,  0.00000000e+00,  5.67765568e-08,
     2.23443223e-05,  2.74725275e-08,  1.39194139e-07,
     5.12820513e-06,  0.00000000e+00,  4.39560440e-08,
    -7.87545788e-06,  0.00000000e+00,  3.47985348e-08],
   [-1.31868132e-03, -7.87545788e-06,  7.87545788e-06,
    -4.17582418e-03, -1.95970696e-05,  5.12820513e-06,
     9.67032967e-03, -2.23443223e-05,  2.23443223e-05,
    -4.17582418e-03, -5.12820513e-06,  1.95970696e-05],
   [ 7.87545788e-06,  3.47985348e-08,  0.00000000e+00,
     1.95970696e-05,  5.67765568e-08,  0.00000000e+00,
    -2.23443223e-05,  1.39194139e-07, -2.74725275e-08,
    -5.12820513e-06,  4.39560440e-08,  0.00000000e+00],
   [-7.87545788e-06,  0.00000000e+00,  3.47985348e-08,
     5.12820513e-06,  0.00000000e+00,  4.39560440e-08,
     2.23443223e-05, -2.74725275e-08,  1.39194139e-07,
    -1.95970696e-05,  0.00000000e+00,  5.67765568e-08],
   [-4.17582418e-03, -1.95970696e-05, -5.12820513e-06,
    -1.31868132e-03, -7.87545788e-06, -7.87545788e-06,
    -4.17582418e-03, -5.12820513e-06, -1.95970696e-05,
     9.67032967e-03, -2.23443223e-05, -2.23443223e-05],
   [ 1.95970696e-05,  5.67765568e-08,  0.00000000e+00,
     7.87545788e-06,  3.47985348e-08,  0.00000000e+00,
    -5.12820513e-06,  4.39560440e-08,  0.00000000e+00,
    -2.23443223e-05,  1.39194139e-07,  2.74725275e-08],
   [-5.12820513e-06,  0.00000000e+00,  4.39560440e-08,
     7.87545788e-06,  0.00000000e+00,  3.47985348e-08,
     1.95970696e-05,  0.00000000e+00,  5.67765568e-08,
    -2.23443223e-05,  2.74725275e-08,  1.39194139e-07]]     
KE1 = np.array(KE1)
xTilde = np.ones(nelx*nely)*VF
xPhys = ((np.tanh(beta*eta)+np.tanh(beta*(xTilde-eta)))/(np.tanh(beta*eta)+np.tanh(beta*(1-eta))))

K = np.zeros((3*(nelx+1)*(nely+1),3*(nelx+1)*(nely+1)))
F =np.zeros((3*(nely+1)*(nelx+1),1))
U1 = np.zeros((3*(nely+1)*(nelx+1),1))
Fe=q*a*b*np.mat([1,b/3,-a/3,1,b/3,a/3,1,-b/3,a/3,1,-b/3,-a/3]).T
xPhys=xPhys.reshape(nelx,nely)

for ely in range(nely):
    for elx in range(nelx):
        n1 = (nely+1)*(elx+1-1)+ely+1
        n2 = (nely+1)*(elx+1) +ely+1
        edof = [3*n1+1, 3*n1+2, 3*n1+3,3*n2+1,3*n2+2,3*n2+3,3*n2-2,3*n2-1, 3*n2, 3*n1-2,3*n1-1, 3*n1]
        edof = [m-1 for m in edof]
        indices = np.ix_(edof, edof)
        K[indices]=K[indices]+(Emin+(1-Emin)*xPhys[ely,elx]**penal)*KE1           
        F[edof,:]=F[edof,:]+Fe   
        
dd = np.where(T2 == 0.01)[0]        
fixeddofs = [(dd)*3+1,(dd)*3+2,(dd)*3+3]
fixeddofs = np.array(fixeddofs)-1  
alldofs = range(3*(nely+1)*(nelx+1))
freedofs =np.setdiff1d(alldofs,fixeddofs).tolist()
K = sparse.csc_matrix(K)
indices = np.ix_(freedofs, freedofs)
U1[freedofs,0]=spsolve(K[indices],F[freedofs,0])
U1[fixeddofs,:]=0

【问题讨论】：

标签： python matlab matrix

【解决方案1】：

假设 K 和 F 的索引对于两个代码都是正确的，差异可能是由于 MATLAB / scipy 各自求解线性方程 Ax = B 的方式。

MATLAB 为其mldivide 函数（又名\）提供了深入的文档。与这种情况最相关的是their algorithm for solving sparse matrices，它可以根据输入矩阵决定使用不同的求解器。在this question thread 中也可以找到关于mldivide 工作原理的很好解释，它将mldivide 与linsolve 进行比较。

同时，scipy 并没有提供太多关于他们在幕后为 spsolve 使用的算法的解释，但似乎他们使用 UMFPACK 作为默认求解器

From scipy's Github, line 189-190:

x = umf.linsolve(umfpack.UMFPACK_A, A, b_vec,autoTranspose=True)

我很难找到 UMFPACK 的文档，但 看来 UMFPACK 然后使用 LU 求解器。 来源：SuiteSparse、Scilab documentation

使用不同的求解器方法将能够解释您的两个代码的数值差异。

因此，根据结果差异的显着程度，您可能需要进行更深入的调查；具体确定您的用例所需的特定求解器，并直接使用这些函数。

【讨论】：

【解决方案2】：

我在这里看到的唯一解释是，在 19684*19684 的情况下，您的 K 矩阵可能几乎是奇异的和/或缩放不良，因此矩阵除法是在最小二乘意义上完成的，或者通过考虑 More-Penrose伪逆或通过最大化解的零分量的数量。 Python 和 Matlab 在这里可能会做出不同的选择。您可以通过计算 K 的数字等级来检查这一点。

【讨论】：