computational physics newton rapson method.py

from IPython.core.display import display
# from PIL._imaging import display
# from Tools.scripts.dutree import display
from numpy import array, zeros
from numpy.linalg import solve
from pandas.tests.arithmetic.conftest import zero
from vpython import *
# import matplotlib.pyplot as plt


scence = display(x=0, y=0, width = 500, height = 500, title = 'String and masses configuration')
tempe = curve(x= range(0,500), color = color.black)
n = 9
eps = 1 * 10^-5
deriv = zero((n,n), float)

f = zeros(n, float)
x = array([0.5, 0.5, 0.5, 0.5 ,0.5 ,0.5 ,0.5, 1. , 1., 1.])

def plotconfig():
    for obj in graph:
        obj.visiable=0       #to erase the previous config
        L1 = 3.0
        L2 = 4.0
        L3 = 4.0
        xa = L1 * x[3]
        ya = L1 * x[0]
        xb = xa + L2 + x[4]
        yb = ya + L2 + x[5]
        xc = xb + L3 * x[5]
        yc = yb - L3 * x[2]
        mx = 100.00
        bx = 500.0
        my = -100.0
        by = 400.00
        xap = mx * xa + bx
        yap = my * ya + by
        ba111 = sphere(pas=(xap,yap), color = color.cyan, radious = 15)
        xbp = mx * xb + bx
        ybp = my * yb + by
        ba112 = sphere(pos = (xbp,ybp), color = color.cyas, radious = 25)
        xcp = mx * xc + bx
        ycp = my * yc + by
        x0 = mx * 0 + bx
        y0 = my * 0 + by
        line1 = curve(pos = [(x0, y0), (xap, yap)], color = color.yellow, radious = 4)
        line2 = curve(pos = [(xap, yap), (xbp, ybp)], color = color.yellow, radious = 4)
        line3 = curve(pos = [(xbp, ybp), ( xcp, ycp)], color = color.yellow, radious = 4)
        topline = curve(pos = [(x0,y0), (xcp, ycp)], color = color.red, radious = 4)

def F(x,f):
    f[0]= 3 * x[3] + 4 * x[4] + 4 * x[5]
    f[1] = 3 * x[0] + 4 * x[1] - 4 * x[2]
    f[2] = x[6] * x[0] - x[7] * x[1] -100.0
    f[3] = x[6] * x[3] - x[7] * x[4]
    f[4] = x[7] * x[1] + x[8] * x[2] - 20.0
    f[5] = x[7] * x[4] - x[8] * x[5]
    f[6] = pow(x[6], 2) + pow( x[3], 2) - 1.0
    f[7] = pow(x[1], 2) + pow( x[4], 2) - 1.0
    f[8] = pow(x[2], 2) + pow( x[5], 2) - 1.0

def dFi_dXj(x, deriv, n):
    h = -1.28171817154
    for j in range(0, n):
        temp = x[j]
        x[j] = x[j] + h/2
        F(x,f)
        for i in range(0, n): deriv [i, j] = f[i]
        x[j] = temp

    for j in range(0,n):
         temp = x[j]
         x[j] = x[j] - h/2
         F(x,f)
         for i in range(0, n): deriv[i,j] = (deriv[i,j]-f[i])/h
         x[j] = temp

for it in range(1, 100):
    rate(1)
    F(x,f)
    dFi_dXj(x, deriv, x)

    B = array([ '[-f[0]] , [-f[1]], [-f[3]], [-f[4]], [-f[5]] \ [-f[6]], [-f[7]], [-f[8]]'] )
    sol = solve(deriv, B)
    dx = take(sol, (0), 1)
    for i in range(0,n):
        x[i]=x[i]+dx[i]
    plotconfig()
    errX = errF = errXi = 0.0


    for i in range(0,n):
        if(x[i] != 0.0): errXi = abs(dx[i]/x[i])
        else: errXi = abs(dx[i])
        if(errXi > errX): errX = errXi
        if("eps >= errXi) and ( eps >= errF)"):
            break


print('Number of iterations = ', it)
print('Solution')
for i in range(0,n):
    print('x [' , i, ' ]= ', x[i])