Сходимость метода дифференциального преобразования (DTM) L/H = 20 ⇐ Python

[Anonymous]

В течение последних нескольких недель я пытался определить естественные частоты вибраций луча Тимошенко в FGM, используя метод дифференциального преобразования (DTM). Чтобы достичь этого, я разработал код Python, где естественные частоты получаются путем решения характерного уравнения с использованием метода Ньютона-Рафсона. Результаты хорошо сходится для соотношения сторон L/H = 5, но сходимость не достигается для l/h = 20.

Код: Выделить всё

import sympy as sp
import numpy as np
from scipy import io, integrate, linalg, signal
from scipy.sparse.linalg import cg, eigs
import sympy as sp
import math
from math import *
from sympy import  *
from sympy import symbols
# Define symbolic variables
d1, d2= symbols('d1 d2 ')
z, x = sp.symbols('z x')
p, Ec, Em, nu = 0, 380., 70., 0.3
rhoc, rhom, alpha = 3960., 2702., 0.
l, Lh, s = 1., 5., 1.
h = l/Lh
b = l / s
ks = 5. / 6.
V1=((1./2.)+(z/h))**p;
Ez1 = (Ec - Em) * V1 + Em - (1/2) * alpha * (Ec + Em)
cnp=integrate(Ez1*z,(z,-h/2,h/2))/(integrate(Ez1,(z,-h/2,h/2)));
print('c=',cnp)
V2 = ((1.  / 2.) + ((z+cnp) / h))**p
# Define functions
Ez = (Ec - Em) * V2 + Em - (1/2) * alpha * (Ec + Em)
rhoz = (rhoc - rhom) *V2 + rhom - (1/2) * alpha * (rhoc + rhom)
fz = z
fz1 = sp.diff(fz, z)
Q11 = Ez
Q55 = Ez / (2 * (1 + nu))
h1=-h/2.-cnp;h2=h/2.-cnp
# Perform symbolic integration

D11 = sp.integrate(Q11 * z**2, (z, h1, h2))
A55a = ks * sp.integrate(Q55 * fz1**2, (z, h1, h2))

I0 = sp.integrate(rhoz, (z, h1, h2))
I2 = sp.integrate(rhoz * z**2, (z, h1, h2))

# Compute coefficients

a11 = (-I0*x) / A55a
a12 = -A55a / A55a
a21 = (-I2*x+ A55a)/D11
a22 = A55a/D11

print("Encastre-Encastre_(CC)")

# Define U and W as symbolic lists
Phi = [0 for _ in range(100)]  # Extend this as needed
W = [0 for _ in range(100)]  # Extend this as needed

# Initialize the given boundary conditions
Phi[0] = 0
Phi[1] = d1
W[0] = 0
W[1] = d2
for N in range(4, 26, 1):
for K in range(0, N):
W[K+2] = (a11 * W[K]+ a12*(K + 1) * Phi[K + 1] ) / ((K+2)*(K+1))
Phi[K+2] = (a21 *  Phi[K] + a22 * (K + 1) * W[K + 1]) /  ((K+2)*(K+1))

# Summation equations
f2 = sum(Phi[k] for k in range(N+1))
f3 = sum(W[k] for k in range(N+1))

f22=sp.expand(f2)
f33=sp.expand(f3)

# Collect terms involving d1, d2
eqt2 = sp.collect(f22, {d1, d2})
eqt3 = sp.collect(f33, {d1, d2})

system = [eqt2, eqt3]

# Extract the coefficients of d1, d2
A = sp.Matrix([[eqt.coeff(d1) for eqt in system],
[eqt.coeff(d2) for eqt in system]])

# Create the vector for the constants
b = sp.Matrix([eqt.subs({d1: 0, d2: 0}) for eqt in system])

# Compute the determinant of the matrix A

det = A.det()
ddet=diff(det,x)
# Convert symbolic expressions to numerical functions
f_numer = sp.lambdify(x, det, 'numpy')
g_numer = sp.lambdify(x, ddet, 'numpy')

# Define functions that return numerical values
def f(x):
return f_numer(x)  # Evaluate determinant at x

def g(x):
return g_numer(x)  # Evaluate derivative at x

# Implementing Newton-Raphson Method
def newtonRaphson(x0, ee, Nst):
print('\n\n*** NEWTON RAPHSON METHOD ***')
step = 1
flag = 1
condition = True

while condition:
if g(x0) == 0.0:  # Avoid division by zero
print('Divide by zero error!')
break

x1 = x0 - f(x0) / g(x0)  # Now properly evaluates numerical values

x0 = x1
step += 1

if step > Nst:
flag = 0
break

condition = abs(f(x1)) > ee  # Ensure numerical comparison

if flag == 1:

print('\nRequired root is: %0.8f' % x1)
Omega_bar=sqrt(x1)*l**2*sqrt(rhom/Em)/h
print(f'N={N:4.2f}   P={p:4.2f}   Lh={Lh:4.2f}   Omega_bar={Omega_bar}')
else:
print('\nNot Convergent.')

# Input Section
x0 =0.009671 #input('Enter Guess: ')
ee = 0.000001#input('Tolerable Error: ')
Nst =1000# input('Maximum Step: ')

# Converting x0 and e to float
x0 = float(x0)
ee = float(ee)

# Converting N to integer
Nst = int(Nst)

#Note: You can combine above three section like this
# x0 = float(input('Enter Guess: '))
# e = float(input('Tolerable Error: '))
# N = int(input('Maximum Step: '))

# Starting Newton Raphson Method

newtonRaphson(x0,ee,Nst)

Пожалуйста, предоставьте любые предложения

Подробнее здесь: https://stackoverflow.com/questions/794 ... dtm-l-h-20