Python与Matlab的矩阵计算差异

Question

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

在最后一步，当我尝试求解位移U= F/K时，出现了一个问题。我已经检查过，对于Matlab和Python，计算出的F和K是相同的。然而，由此产生的U是different.When矩阵大小达到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

Answer 1

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

MATLAB为他们的mldivide函数提供了深入的文档(也称为.\)。与这种情况最相关的是their algorithm for solving sparse matrices，它可以根据输入矩阵决定使用不同的求解器。关于mldivide如何工作的一个很好的解释也可以在this question thread中找到，它将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

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

因此，根据结果中的差异有多大，您可能希望进行更深入的研究；具体决定您的用例需要哪些特定的求解器，并直接使用这些函数。

Answer 2

我在这里看到的唯一解释是，在19684*19684的情况下，您的K矩阵可能几乎是奇异的和/或缩放得很差，因此矩阵除法是在最小二乘意义上完成的，要么考虑更多的彭罗斯伪逆，要么通过最大化解决方案的零分量的数量。Python和Matlab在这里可能会做出不同的选择。你可以通过计算K的数值秩来检查这一点。