[已解决]关于sympy和numpy的使用问题：can't convert expression to float

xxfholic

问题描述：本人Python初学者，最近在完成老师所写的大作业，写了一段代码，发现在使用numpy数组时在  M[]  处出现can't convert expression to float的错误，希望大佬斧正，如果顺便帮我检查有没有其他错误就更好了，跪谢！

传递矩阵法
import numpy as np
import sympy as sym

E = 2.06e11
I = np.array([
    0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08
])
d = np.array([
    0.025, 0.025, 0.025, 0.025, 0.025, 0.025, 0.025, 0.025, 0.025
])
m = np.array([
    0.308, 0.308, 9.21, 0.308, 0.308, 0.308, 0.308, 0.308, 0
])
Id = np.array([
    0, 0, 0.03196, 0, 0, 0, 0, 0, 0
])
Io = np.array([
    0, 0, 0.06392, 0, 0, 0, 0, 0, 0
])
J = np.array([
    1.917e-8 ,1.917e-8 ,1.917e-8 ,1.917e-8 ,1.917e-8 ,1.917e-8 ,1.917e-8 ,1.917e-8 ,1.917e-8
])


L = np.zeros([9,4,4])
for i in np.arange(0,9):
    L[i] = np.array([
        [ 1 , I[i], I[i]**2/2/E/J[i], I[i]**3/6/E/J[i] ],
        [ 0, 1, I[i]/E/J[i], I[i]**2/2/E/J[i] ],
        [ 0, 0, 1, I[i] ],
        [ 0, 0, 0, 1 ]
    ])

w = sym.symbols('w')

M = np.zeros([9,4,4])
for i in np.arange(0,9):
    if i <= 1 or i >= 3:
        M[i] = np.array([
            [ 1, 0, 0, 0 ],
            [ 0, 1, 1, 0 ],
            [ 0, 0, 1, 0 ],
            [ m[i]*w**2, 0, 0, 1 ]
        ])
    else:
        M[i] = np.array([
            [ 1, 0, 0, 0],
            [ 0, 1, 0, 0 ],
            [0, (Io[i]/Id[i]-1)*Id[i]*w**2, 1, 0 ],
            [ m[i]*w**2, 0, 0, 1 ]
        ])

#print(M)
H = np.zero([4,4])
for i in np.arange(0,9):
    H[i] *= M[i] * L[i]

h12, h14, h34, h32 = H[1,2], H[1,4], H[3,4], H[3,2]
delta = h12 * h34 - h14 * h32

sym.solve(delta-w,w)
print(w)

最佳答案

月排行榜 / 总排行榜

永恒的蓝色梦想

2020-4-17 18:32:36

m = np.array([
    0.308, 0.308, 9.21, 0.308, 0.308, 0.308, 0.308, 0.308, 0
])

w = sym.symbols('w')

M[i] = np.array([
            [ 1, 0, 0, 0 ],
            [ 0, 1, 1, 0 ],
            [ 0, 0, 1, 0 ],
            [ m[i]*w**2, 0, 0, 1 ]
        ])

跨库调用，有点NB

不是什么库都兼容的啊

跳转到最佳答案楼层

xxfholic

附上原来的MATLAB代码

%Rotor critical speed calculation 转子临界转速计算
%*****************************************************************************
%转子数据
E=2.06e11;%轴的弹性模量 N/m2
l=[0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08];%轴总长为 72cm,分成 9 段每段的长度相等都是 0.08m
d=[0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025];%轴的直径 m
m=[0.308,0.308,9.21,0.308,0.308,0.308,0.308,0.308,0];%各段的质量
Id=[0,0,0.03196,0,0,0,0,0,0];% 直 径 转 动 惯 量 kg*m2 
Io=[0,0,0.06392,0,0,0,0,0,0];%极转动惯量
J=[1.917e-8,1.917e-8,1.917e-8,1.917e-8,1.917e-8,1.917e-8,1.917e-8,1.917e-8,1.917e-8];%J 为各段面积惯性矩，J=pi*d^4/64
% E=2.06e7 %轴的弹性模量 N/cm2
% l=[8,8,8,8,8,8,8,8,8];%轴总长为 72cm,分成 9 段每段的长度相等都是 8cm
% d=[2.5,2.5,2.5,2.5,2.5,2.5,2.5,2.5,2.5];%轴的直径
% m=[0.308,0.308,9.21,0.308,0.308,0.308,0.308,0.308,0];%各段的质量
% Id=[0,0,319.6,0,0,0,0,0,0];%直径转动惯量
% Io=[0,0,639.2,0,0,0,0,0,0];%极转动惯量
% J=[1.917,1.917,1.917,1.917,1.917,1.917,1.917,1.917,1.917];%J 为各段面积惯性矩，J=pi*d^4/64
%*****************************************************************************
%1-9 段的传递矩阵 L clc
L=cell(1,9);
for i=1:9
L{1,i}=[1,l(i),(l(i))^2/(2*E*J(i)),(l(i))^3/(6*E*J(i));
0,        1,        l(i)/(E*J(i)),        (l(i))^2/(2*E*J(i));
0,        0,        1,        l(i);
0,        0,        0,        1];
end
%1-9 站的传递矩阵 M M=cell(1,9);
syms w;
 
for i=1:9

if i<=2||i>=4
M{1,i}=[1,        0,        0,        0;
0,        1,        0,        0;
0,        0,        1,        0;
m(i)*w^2,        0,        0,        1];
else
M{1,i}=[1,        0,        0,        0;
0,        1,        0,        0;
0, ((Io(i)/Id(i))-1)*Id(i)*w^2, 1,        0;
m(i)*w^2,        0,        0,        1];
end
end
%*****************************************************************************
%传递矩阵计算
%左右两端的边界条件,转子两边看成是是简支的
% BJ9=[0,Seta9,0,Q9];
% BJ0=[0,Seta0,0,Q0];
% BJ9=M{1,9}*L{1,9}*M{1,8}*L{1,8}*M{1,7}*L{1,7}*M{1,6}*L{1,6}*M{1,5}*L{1,5}
% *M{1,4}*L{1,4}*M{1,3}*L{1,3}*M{1,2}*L{1,2}*M{1,1}*L{1,1}*BJ0
H4x4=M{1,9}*L{1,9}*M{1,8}*L{1,8}*M{1,7}*L{1,7}*M{1,6}*L{1,6}*M{1,5}*L{1,5}*...
M{1,4}*L{1,4}*M{1,3}*L{1,3}*M{1,2}*L{1,2}*M{1,1}*L{1,1}; h12=H4x4(1,2);h14=H4x4(1,4);h32=H4x4(3,2);h34=H4x4(3,4);
%采用作图法，根据传递矩阵法的原理，画出Δ随Ω（轴的转速）
%变化的曲线，曲线与Δ=0 时的水平轴交点频率就为临界转速。
Drta=h12*h34-h14*h32;
% w=200:1:1600;
% plot(w,subs(Drta,w))%作图观察第一阶临界转速的变化趋势和 Drta=0 轴的交点
%        xlabel('转速/rad/s','FontSize',10.5);ylabel('Δ值','FontSize',10.5);
% grid on
%*****************************************************************************
% a=subs(Drta,w)
% Drta=vpa(a)
 
%采用试凑法，选定一系列的w 值，代入 Drta 式子中进行计算，得到一系列的 Drta 值，当Drta 值满足一定的
%的精度时就可以认为此时的 W 值为临界转数。
for w=200:1:400
a=subs(Drta,w);
if        a>=-0.01&&a<=0.001%计算精度的选取很重要，精度不能太小，否则得不到准确结果w1cor=w
break
end
end
%注意利用传递矩阵发计算转子的临界速度对 MATLAB 的计算精度有一定的要求，所以在计
%算中一定要采用双精度数，且各个计算物理量都要采用国际单位,才能得到正确的临界转数
%*****************************************************************************

xxfholic

xxfholic

永恒的蓝色梦想

m = np.array([
    0.308, 0.308, 9.21, 0.308, 0.308, 0.308, 0.308, 0.308, 0
])

w = sym.symbols('w')

M[i] = np.array([
            [ 1, 0, 0, 0 ],
            [ 0, 1, 1, 0 ],
            [ 0, 0, 1, 0 ],
            [ m[i]*w**2, 0, 0, 1 ]
        ])

跨库调用，有点NB

不是什么库都兼容的啊

xxfholic

永恒的蓝色梦想 发表于 2020-4-17 18:32
跨库调用，有点NB

不是什么库都兼容的啊

我也是想过这个问题，可是如果用Python的话，还能有别的办法去解决这个问题吗？谢谢大佬指点！

永恒的蓝色梦想

xxfholic 发表于 2020-4-17 19:39
我也是想过这个问题，可是如果用Python的话，还能有别的办法去解决这个问题吗？谢谢大佬指点！

你可倒是把题目放出来啊