Python Forum
how to solve complex equations in python
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
how to solve complex equations in python
#1
Heart 
I am trying to find the solution to the equation. Det(matrix A)=0, var:omega
###################################################################

from sympy import  Function, symbols,atan, sin, cos, pi, Matrix, solve

A11, A12, A22, A66, B11, B12, B22, B66, C11, C12, C22, C66, D11, D12, D22, D66, E11, E12, E22, E66, G11, G12, G22, G66, F44, F55, I1, I2, I3, I4, I5, I6\
   =symbols('A11 A12 A22 A66 B11 B12 B22 B66 C11 C12 C22 C66 D11 D12 D22 D66 E11 E12 E22 E66 G11 G12 G22 G66 F44 F55 I1 I2 I3 I4 I5 I6')    

omega =symbols('omega')

# from sympy.abc import x, y, t

u, U, v, V, w, W ,phix, Phix, phiy, Phiy, E, rho, muy, fz, gz= symbols('u U v V w W phix Phix phiy Phiy E rho muy fz gz', cls=Function)

from scipy.integrate import quad

###############################################################################
h=0.002  
a=20*h
b=a
#mode
m=1
n=1

k0=0
k1=0
###############################################################################
def fz(z):
    return h*(atan(3*z/h)+sin(2*pi*z/h)/3/pi)/5

def gz(z):
    return (1/(1+2*(z**3)/(h**3))+cos(3*pi*z/h))/5
###############################################################################
Em=70*(10**9)
Ec=380*(10**9)
rhom=2702
rhoc=3800
muym=0.3
muyc=0.3

alpha=0
p=1
###############################################################################
def E(z):
    return (Ec-Em)*((0.5+z/h)**p)+Em-alpha*(Ec+Em)/2
def rho(z):
    return (rhoc-rhom)*((0.5+z/h)**p)+rhom-alpha*(rhoc+rhom)/2
def muy(z):
    return (muyc-muym)*((0.5+z/h)**p)+muym-alpha*(muyc+muym)/2
###############################################################################
#Ma tran do cung Q
def Q11(z):
    return E(z)/(1-muy(z)**2)
def Q12(z):
    return E(z)*muy(z)/(1-muy(z)**2)
def Q22(z):
    return E(z)/(1-muy(z)**2)
def Q66(z):
    return E(z)/(1+muy(z))/2
###############################################################################
A11= quad(lambda z: Q11(z), -h/2, h/2)[0]
A12= quad(lambda z: Q12(z), -h/2, h/2)[0]
A22= quad(lambda z: Q22(z), -h/2, h/2)[0]
A66= quad(lambda z: Q66(z), -h/2, h/2)[0]

B11= quad(lambda z: Q11(z)*z, -h/2, h/2)[0]
B12= quad(lambda z: Q12(z)*z, -h/2, h/2)[0]
B22= quad(lambda z: Q22(z)*z, -h/2, h/2)[0]
B66= quad(lambda z: Q66(z)*z, -h/2, h/2)[0]

C11= quad(lambda z: Q11(z)*fz(z), -h/2, h/2)[0]
C12= quad(lambda z: Q12(z)*fz(z), -h/2, h/2)[0]
C22= quad(lambda z: Q22(z)*fz(z), -h/2, h/2)[0]
C66= quad(lambda z: Q66(z)*fz(z), -h/2, h/2)[0]

D11= quad(lambda z: Q11(z)*z*z, -h/2, h/2)[0]
D12= quad(lambda z: Q12(z)*z*z, -h/2, h/2)[0]
D22= quad(lambda z: Q22(z)*z*z, -h/2, h/2)[0]
D66= quad(lambda z: Q66(z)*z*z, -h/2, h/2)[0]

E11= quad(lambda z: Q11(z)*fz(z)*z, -h/2, h/2)[0]
E12= quad(lambda z: Q12(z)*fz(z)*z, -h/2, h/2)[0]
E22= quad(lambda z: Q22(z)*fz(z)*z, -h/2, h/2)[0]
E66= quad(lambda z: Q66(z)*fz(z)*z, -h/2, h/2)[0]

G11= quad(lambda z: Q11(z)*fz(z)*fz(z), -h/2, h/2)[0]
G12= quad(lambda z: Q12(z)*fz(z)*fz(z), -h/2, h/2)[0]
G22= quad(lambda z: Q22(z)*fz(z)*fz(z), -h/2, h/2)[0]
G66= quad(lambda z: Q66(z)*fz(z)*fz(z), -h/2, h/2)[0]

F44= quad(lambda z: Q66(z)*gz(z)*gz(z), -h/2, h/2)[0]
F55= quad(lambda z: Q66(z)*gz(z)*gz(z), -h/2, h/2)[0]

I1=quad(lambda z: rho(z), -h/2, h/2)[0]
I2=quad(lambda z: rho(z)*z, -h/2, h/2)[0]
I3=quad(lambda z: rho(z)*z*z, -h/2, h/2)[0]
I4=quad(lambda z: rho(z)*fz(z), -h/2, h/2)[0]
I5=quad(lambda z: rho(z)*fz(z)*z, -h/2, h/2)[0]
I6=quad(lambda z: rho(z)*fz(z)*fz(z), -h/2, h/2)[0]




a11=(a*A11*b**2*m**2*pi**2+a**3*A66*n**2*pi**2)/(a**3*b**2);
a12=(a**2*A12*b*m*n*pi**2+a**2*A66*b*m*n*pi**2)/(a**3*b**2);
a13=((-b**2)*B11*m**3*pi**3-a**2*B12*m*n**2*pi**3-2*a**2*B66*m*n**2*pi**3)/(a**3*b**2);
a14=((b**2*C11*m**2+a**2*C66*n**2)*pi**2)/(a**2*b**2);
a15=((C12+C66)*m*n*pi**2)/(a*b);
a21=(a*A12*b**2*m*n*pi**2+a*A66*b**2*m*n*pi**2)/(a**2*b**3);
a22=(A66*b**3*m**2*pi**2+a**2*A22*b*n**2*pi**2)/(a**2*b**3);
a23=((-b**2)*B12*m**2*n*pi**3-2*b**2*B66*m**2*n*pi**3-a**2*B22*n**3*pi**3)/(a**2*b**3);
a24=((C12+C66)*m*n*pi**2)/(a*b);
a25=((b**2*C66*m**2+a**2*C22*n**2)*pi**2)/(a**2*b**2);

a31=((-a)*b**4*B11*m**3*pi**3-a**3*b**2*B12*m*n**2*pi**3-2*a**3*b**2*B66*m*n**2*pi**3)/(a**4*b**4);
a32=((-a**2)*b**3*B12*m**2*n*pi**3-2*a**2*b**3*B66*m**2*n*pi**3-a**4*b*B22*n**3*pi**3)/(a**4*b**4);
a33=(1/(a**4*b**4))*(a**4*b**4*k0+a**2*b**4*k1*m**2*pi**2+a**4*b**2*k1*n**2*pi**2+b**4*D11*m**4*pi**4+2*a**2*b**2*D12*m**2*n**2*pi**4+4*a**2*b**2*D66*m**2*n**2*pi**4+a**4*D22*n**4*pi**4);
a34=-((m*(b**2*E11*m**2+a**2*(E12+2*E66)*n**2)*pi**3)/(a**3*b**2));
a35=-((n*(b**2*(E12+2*E66)*m**2+a**2*E22*n**2)*pi**3)/(a**2*b**3));

a41=(a*b**2*C11*m**2*pi**2+a**3*C66*n**2*pi**2)/(a**3*b**2);
a42=(a**2*b*C12*m*n*pi**2+a**2*b*C66*m*n*pi**2)/(a**3*b**2);
a43=((-b**2)*E11*m**3*pi**3-a**2*E12*m*n**2*pi**3-2*a**2*E66*m*n**2*pi**3)/(a**3*b**2);
a44=(a*b**2*G11*m**2*pi**2+a**3*(b**2*F55+G66*n**2*pi**2))/(a**3*b**2);
a45=((G12+G66)*m*n*pi**2)/(a*b);

a51=(a*b**2*C12*m*n*pi**2+a*b**2*C66*m*n*pi**2)/(a**2*b**3);
a52=(b**3*C66*m**2*pi**2+a**2*b*C22*n**2*pi**2)/(a**2*b**3);
a53=((-b**2)*E12*m**2*n*pi**3-2*b**2*E66*m**2*n*pi**3-a**2*E22*n**3*pi**3)/(a**2*b**3);
a54=((G12+G66)*m*n*pi**2)/(a*b);
a55=(b**3*G66*m**2*pi**2+a**2*(b**3*F44+b*G22*n**2*pi**2))/(a**2*b**3);



A=Matrix([[a11-I1*omega,a12,a13+(I2*m*pi*omega)/a,a14-I4*omega,a15],\
          [a21,a22-I1*omega,a23+(I2*n*pi*omega)/b,a24,a25-I4*omega],\
          [a31+(I2*m*pi*omega)/a,a32+(I2*n*pi*omega)/b,a33-((a**4*b**4*I1+a**2*b**4*I3*m**2*pi**2+a**4*b**2*I3*n**2*pi**2)*omega)/(a**4*b**4),a34+(I5*m*pi*omega)/a,a35+(I5*n*pi*omega)/b],\
          [a41-I4*omega,a42,a43+(I5*m*pi*omega)/a,a44-I6*omega,a45],\
          [a51,a52-I4*omega,a53+(I5*n*pi*omega)/b,a54,a55-I6*omega]])
    
DET=A.det()

ng=solve(DET, omega)

print(ng)
Larz60+ write Dec-26-2021, 09:44 AM:
Please post all code, output and errors (it it's entirety) between their respective tags. Refer to BBCode help topic on how to post. Use the "Preview Post" button to make sure the code is presented as you expect before hitting the "Post Reply/Thread" button.
Fixed for you this time. Please use bbcode tags on future posts.
Reply
#2
Help me with code to solve complex equations in python........... Thank you so much <33333
Reply
#3
What are you trying to do?
Reply
#4
(Dec-26-2021, 02:18 PM)Linenloid Wrote: What are you trying to do?
I am trying to find the solution to the equation. Det(matrix A)=0, var:omega. Thank you so much
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
  Python selenium + Xevil to solve recaptcha sunny9495 6 3,284 Apr-04-2022, 10:23 AM
Last Post: sunny9495
  SOlving LInear Equations in Python(Symoy, NUmpy) - Coefficient Problem quest 3 1,689 Jan-30-2022, 10:53 PM
Last Post: quest
Question Help with code to solve polynomials equations hiviera 1 1,761 Jul-31-2021, 01:56 AM
Last Post: hiviera
  NotImplementedError: pseudo-class is not implemented - how to Update Python to solve apollo 1 3,026 May-16-2021, 08:03 AM
Last Post: buran
  Solve simple equation in Python kmll 4 3,024 Nov-01-2020, 04:34 PM
Last Post: deanhystad
  Solve system of equations Sancho_Pansa 19 8,692 Oct-27-2020, 08:15 AM
Last Post: Sancho_Pansa
  How to solve equations, with groups of variables and or constraints? ThemePark 0 1,644 Oct-05-2020, 07:22 PM
Last Post: ThemePark
  Solve a system of linear equations with binary variables lopeslimagabriel 3 2,443 Sep-24-2020, 07:09 AM
Last Post: scidam
  Differential equations with initial condition in Python (change a working code) Euler2 1 1,793 May-29-2020, 04:06 PM
Last Post: Euler2
  How to solve difficult non-linear equations? shreeniket 3 2,348 Apr-23-2020, 01:36 AM
Last Post: shreeniket

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020