[Python] Convergence of differential transformation method (DTM) L/h=20

For the past few weeks, I have been trying to determine the natural frequencies of the vibrations of a Timoshenko beam in FGM using the differential transformation method (DTM). To achieve this, I developed a Python code where the natural frequencies are obtained by solving a characteristic equation using the Newton-Raphson method. The results converge well for an aspect ratio of L/h=5, but convergence is not achieved for 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)


Please provide any suggestions

Continue reading...