Python中文网 - 问答频道, 解决您学习工作中的Python难题和Bug
Python常见问题
我正在将一个有限元分析代码从Matlab传输到Python
当我试图求解位移U=F/K时，在最后一步出现了一个问题。 我已经检查了计算出的F和K对于Matlab和Python是相同的。但是结果U是不同的。当矩阵大小达到19684*19684时，有效数字将完全不同，甚至符号将相反
我的问题是，对于相同的K和F矩阵，为什么Python和MATLAB会对U=F/K给出不同的答案。需要注意的是，如果K和F的大小很小（例如1000*10001000*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
Tags： fornpzerosnubetaetaelyn2
热门问题
- 如何添加虚拟方法
- 如何添加表示整数的擦边字符串？
- 如何添加要在Bokeh中使用的新font.ttf文件？
- 如何添加要显示的矩阵XY轴编号和XY轴
- 如何添加计数？
- 如何添加计数器函数？
- 如何添加计数器列来计算数据帧中另一列中的特定值？
- 如何添加计数器来跟踪while循环中的月份和年份？
- 如何添加计数并删除countplot的顶部和右侧脊椎？
- 如何添加计时器wx.应用程序更新窗口对象的主循环？
- 如何添加评论到帖子？PostDetailVew，Django 2.1.5
- 如何添加评论拉梅尔亚姆
- 如何添加诸如矩阵Python/Pandas之类的数据帧？
- 如何添加谷歌地点自动完成到Flask？
- 如何添加超时、python discord bot
- 如何添加超过1dp的检查
- 如何添加距离方法
- 如何添加跟随游戏的敌人精灵
- 如何添加路径以便python可以找到程序？
- 如何添加身份验证/安全性以使用happybase访问HBase？
热门文章
- Python覆盖写入文件
- 怎样创建一个 Python 列表？
- Python3 List append()方法使用
- 派森语言
- Python List pop()方法
- Python Django Web典型模块开发实战
- Python input() 函数
- Python3 列表(list) clear()方法
- Python游戏编程入门
- 如何创建一个空的set?
- python如何定义(创建)一个字符串
- Python标准库 [The Python Standard Library by Ex
- Python网络数据爬取及分析从入门到精通（分析篇）
- Python3 for 循环语句
- Python List insert() 方法
- Python 字典(Dictionary) update()方法
- Python编程无师自通 专业程序员的养成
- Python3 List count()方法
- Python 网络爬虫实战 [Web Crawler With Python]
- Python Cookbook（第2版）中文版
假设
K和F的索引对于这两种代码都是正确的，那么差异可能是由于MATLAB/scipy各自如何解决线性方程Ax = BMATLAB为它们的^{} 函数（又名} 进行了比较
\）提供了深入的文档。与这种情况最相关的是their algorithm for solving sparse matrices，它可以根据输入矩阵决定使用不同的解算器。关于mldivide如何工作的一个很好的解释也可以在this question thread中找到，它将mldivide与^{同时，scipy并没有提供太多关于他们在幕后对^{} 使用什么算法的解释，但似乎他们使用UMFPACK作为默认解算器
From scipy's Github, line 189-190:
我很难找到UMFPACK的文档，但是似乎UMFPACK然后使用LU解算器来源：SuiteSparse，Scilab documentation
使用不同的解算器方法将能够解释两个代码的数值差异
因此，根据结果差异的显著程度，您可能需要进行更深入的调查；具体决定您的用例需要什么特定的解算器，并直接使用这些函数
我在这里看到的唯一解释是，在19684*19684的情况下，你的K矩阵可能几乎是奇异的和/或比例很差，因此矩阵除法是在最小二乘意义下进行的，或者通过考虑Penrose伪逆，或者通过最大化解的零分量数。Python和Matlab在这里可能会做出不同的选择。您可以通过计算K的数值秩来检查这一点
相关问题 更多 >
编程相关推荐