Код: Выделить всё
import sympy as sp
import numpy as np
from scipy import io, integrate, linalg, signal
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
from sympy import lambdify
# Define symbolic variable
# Define symbols for c1, c2, c3, and other variables
x,c1, c2, c3, k = symbols('x c1 c2 c3 k')
A11=4.9091
B11=0.0294
D11=0.0014
I0=140.8182
I1=0.1191
I2=0.0307
# Define expressions for e1, e2, e3, e4, e5, and e6
e1 = -I0 * x / A11
e2 = I1 * x / A11
e3 = B11 / A11
e4 = -(B11 * e1 + I1 * x) / (B11 * e3 - D11)
e5 = -I0 * x / (B11 * e3 - D11)
e6 = (-B11 * e2 + I2 * x) / (B11 * e3 - D11)
# Define U and W as symbolic lists
U = [0 for _ in range(50)] # Extend this as needed
W = [0 for _ in range(50)] # Extend this as needed
# Initialize the given boundary conditions
U[0] = 0
U[1] = c1
W[0] = 0
W[1] = 0
W[2] = c2
W[3] = c3
# Loop over N and K, as in the Maple code
for N in range(12, 30, 2):
for K in range(0, N):
U[K+2] = (e1*U[K] + e2*(K+1)*W[K+1] + e3*(K+3)*(K+2)*(K+1)*W[K+3]) / ((K+2)*(K+1))
W[K+4] = (e4*(K+1)*U[K+1] + e5*W[K] + e6*(K+2)*(K+1)*W[K+2]) / ((K+4)*(K+3)*(K+2)*(K+1))
# Define the sum equations
f1 = sum(U[k] for k in range(N+1)) # Sum for U
f2 = sum(W[k] for k in range(N+1)) # Sum for W
f3 = sum(k * W[k] for k in range(N+1)) # Sum for W weighted by k
f11=sp.expand(f1)
f22=sp.expand(f2)
f33=sp.expand(f3)
# Collect terms involving c1, c2, c3
eqt1 = sp.collect(f11, {c1, c2, c3})
eqt2 = sp.collect(f22, {c1, c2, c3})
eqt3 = sp.collect(f33, {c1, c2, c3})
system = [eqt1, eqt2, eqt3]
# Extract the coefficients of c1, c2, c3
A = sp.Matrix([[eqt.coeff(c1) for eqt in system],
[eqt.coeff(c2) for eqt in system],
[eqt.coeff(c3) for eqt in system]])
# Create the vector for the constants
b = sp.Matrix([eqt.subs({c1: 0, c2: 0, c3: 0}) for eqt in system])
# Compute the determinant of the matrix A
det = A.det()
ddet=diff(det,x)
# Defining Function
def f(x):
return det
# Defining derivative of function
def g(x):
return ddet
# Implementing Newton Raphson Method
def newtonRaphson(x0,e,Nst):
print('\n\n*** NEWTON RAPHSON METHOD IMPLEMENTATION ***')
step = 1
flag = 1
condition = True
while condition:
if g(x0) == 0.0:
print('Divide by zero error!')
break
x1 = x0 - f(x0)/g(x0)
#print('Iteration-%d, x1 = %0.6f and f(x1) = %0.6f' % (step, x1, f(x1)))
x0 = x1
step = step + 1
if step > Nst:
flag = 0
break
condition = abs(f(x1)) > e
if flag==1:
print('\nRequired root is: %0.8f' % x1)
else:
print('\nNot Convergent.')
# Input Section
x0 =0.1 #input('Enter Guess: ')
e = 0.0001#input('Tolerable Error: ')
Nst =10# input('Maximum Step: ')
# Converting x0 and e to float
x0 = float(x0)
e = float(e)
# 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,e,Nst)
Подробнее здесь: https://stackoverflow.com/questions/794 ... on-raphson