notebook.community

Edit and run



In [1]:

    
# import numpy as np

# # !/usr/bin/env python3
# # -*- coding: utf-8 -*-
# """
# Created on 20181219

# @author: zhangji

# Trajection of a ellipse, Jeffery equation. 
# """

# %pylab inline
# pylab.rcParams['figure.figsize'] = (25, 11)
# fontsize = 40

# import numpy as np
# import scipy as sp
# from scipy.optimize import leastsq, curve_fit
# from scipy import interpolate
# from scipy.interpolate import interp1d
# from scipy.io import loadmat, savemat
# # import scipy.misc

# import matplotlib
# from matplotlib import pyplot as plt
# from matplotlib import animation, rc
# import matplotlib.ticker as mtick
# from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes
# from mpl_toolkits.mplot3d import Axes3D, axes3d

# from sympy import symbols, simplify, series, exp
# from sympy.matrices import Matrix
# from sympy.solvers import solve

# from IPython.display import display, HTML
# from tqdm import tqdm_notebook as tqdm
# import pandas as pd
# import re
# from scanf import scanf
# import os
# import glob

# from codeStore import support_fun as spf
# from src.support_class import *
# from src import stokes_flow as sf

# rc('animation', html='html5')
# PWD = os.getcwd()
# font = {'size': 20}
# matplotlib.rc('font', **font)
# np.set_printoptions(linewidth=90, precision=5)

import numpy as np
from sympy import *
from sympy import lambdify
from sympy.utilities.lambdify import lambdastr

import matplotlib
from matplotlib import pyplot as plt
from matplotlib import animation, rc
import matplotlib.ticker as mtick
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes
from mpl_toolkits.mplot3d import Axes3D, axes3d

from time import time
from src.support_class import *
from src import jeffery_model as jm
from IPython.display import display

init_printing(use_latex='mathjax')



In [5]:

    
x0, x1, x2 = symbols('x0:3')
lbd = Symbol('lbd')

# u0 = Function('u0')(x0, x1, x2)
# u1 = Function('u1')(x0, x1, x2)
# u2 = Function('u2')(x0, x1, x2)

# u0 = Function('u0')(x1, x2)
# u1 = Function('u1')(x2, x0)
# u2 = 0

# f0, f1, f2 = symbols('f0:3')
# X = Matrix([x0, x1, x2])
# r = ((X.T * X) ** 0.5) [0]
# F = Matrix([f0, f1, f2])
# # F = Matrix([1, 0, 0])
# Mij = eye(3) / r + X * X.T / (r ** 3)
# u0, u1, u2 = simplify(Mij * F)

f0, f1, f2 = symbols('f0:3')
h = symbols('h')
X = Matrix([x0, x1, x2])
dX1 = Matrix([x0, x1, x2 - h])
dX2 = Matrix([x0, x1, x2 + h])
r1 = ((dX1.T * dX1) ** 0.5) [0]
r2 = ((dX2.T * dX2) ** 0.5) [0]
F = Matrix([f0, f1, f2])
Q = Matrix([-f0, -f1, f2])
B = Matrix([-f0, -f1, f2]).T
G = Matrix([0, 0, -1])
M1ij = 1 / (8 * pi) * (eye(3) / r1 + dX1 * dX1.T / (r1 ** 3))
M2ij = 1 / (8 * pi) * (eye(3) / r2 + dX2 * dX2.T / (r2 ** 3))
Mpdij = 1 / (4 * pi) * (eye(3) / r2 ** 3 - 3 * dX2 * dX2.T / (r2 ** 5))
us10, us11, us12 = simplify(M1ij * F)
us20, us21, us22 = simplify(M2ij * F)
upd0, upd1, upd2 = simplify(Mpdij * Q)
tu0, tu1, tu2 = simplify(M2ij * G)
usd0, usd1, usd2 = simplify(B * Matrix([[diff(tu0, x0), diff(tu0, x1), diff(tu0, x2)], 
                                        [diff(tu1, x0), diff(tu1, x1), diff(tu1, x2)], 
                                        [diff(tu2, x0), diff(tu2, x1), diff(tu2, x2)]]).T)
u0 = simplify(us10 - us20 - 2 * h * usd0 - h ** 2 * upd0)
u1 = simplify(us11 - us21 - 2 * h * usd1 - h ** 2 * upd1)
u2 = simplify(us12 - us22 - 2 * h * usd2 - h ** 2 * upd2)
# Jij = simplify(Matrix([[diff(u0, x0), diff(u0, x1), diff(u0, x2)], 
#               [diff(u1, x0), diff(u1, x1), diff(u1, x2)], 
#               [diff(u2, x0), diff(u2, x1), diff(u2, x2)]]))
# Sij = simplify(1 / 2 * (Jij + Jij.T))
# Oij = simplify(1 / 2 * (Jij - Jij.T))
display(u0)

# print('Trace B')
# Bij = simplify(Oij + lbd * Sij)
# display(simplify(Bij[0, 0] + Bij[1, 1] + Bij[2, 2]))

# print()
# print('Trace B2')
# Bij2 = simplify(Bij*Bij)
# display(simplify(Bij2[0, 0] + Bij2[1, 1] + Bij2[2, 2]))

# print()
# print('Trace B3')
# Bij3 = simplify(Bij*Bij*Bij)
# display(Bij3[0, 0] + Bij3[1, 1] + Bij3[2, 2])









    





$$\frac{0.75 f_{0} h x_{0}^{2} x_{2}}{\pi \left(x_{0}^{2} + x_{1}^{2} + \left(h + x_{2}\right)^{2}\right)^{2.5}} - \frac{f_{0} h x_{2}}{4 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h + x_{2}\right)^{2}\right)^{1.5}} + \frac{f_{0} x_{0}^{2}}{8 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h - x_{2}\right)^{2}\right)^{1.5}} - \frac{f_{0} x_{0}^{2}}{8 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h + x_{2}\right)^{2}\right)^{1.5}} + \frac{f_{0}}{8 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h - x_{2}\right)^{2}\right)^{0.5}} - \frac{f_{0}}{8 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h + x_{2}\right)^{2}\right)^{0.5}} + \frac{0.75 f_{1} h x_{0} x_{1} x_{2}}{\pi \left(x_{0}^{2} + x_{1}^{2} + \left(h + x_{2}\right)^{2}\right)^{2.5}} + \frac{f_{1} x_{0} x_{1}}{8 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h - x_{2}\right)^{2}\right)^{1.5}} - \frac{f_{1} x_{0} x_{1}}{8 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h + x_{2}\right)^{2}\right)^{1.5}} + \frac{3 f_{2} h^{3} x_{0}}{4 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h + x_{2}\right)^{2}\right)^{2.5}} + \frac{3 f_{2} h^{2} x_{0} x_{2}}{4 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h + x_{2}\right)^{2}\right)^{2.5}} - \frac{0.75 f_{2} h x_{0} \left(h + x_{2}\right)^{2}}{\pi \left(x_{0}^{2} + x_{1}^{2} + \left(h + x_{2}\right)^{2}\right)^{2.5}} - \frac{f_{2} h x_{0}}{8 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h - x_{2}\right)^{2}\right)^{1.5}} + \frac{f_{2} h x_{0}}{8 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h + x_{2}\right)^{2}\right)^{1.5}} + \frac{f_{2} x_{0} x_{2}}{8 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h - x_{2}\right)^{2}\right)^{1.5}} - \frac{f_{2} x_{0} x_{2}}{8 \pi \left(x_{0}^{2} + x_{1}^{2} + \left(h + x_{2}\right)^{2}\right)^{1.5}}$$



In [31]:

    
# lam_u0 = lambdify((x0, x1, x2, f0, f1, f2, h), u0)
print(lambdastr((x0, x1, x2, f0, f1, f2, h), u0))
print()
print(lambdastr((x0, x1, x2, f0, f1, f2, h), u1))
print()
print(lambdastr((x0, x1, x2, f0, f1, f2, h), u2))
print()









    



lambda x0,x1,x2,f0,f1,f2,h: (0.75*f0*h*x0**2*x2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi - 1/4*f0*h*x2*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + (1/8)*f0*x0**2*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi - 1/8*f0*x0**2*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + (1/8)*f0*(x0**2 + x1**2 + (h - x2)**2)**(-0.5)/math.pi - 1/8*f0*(x0**2 + x1**2 + (h + x2)**2)**(-0.5)/math.pi + 0.75*f1*h*x0*x1*x2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi + (1/8)*f1*x0*x1*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi - 1/8*f1*x0*x1*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + (3/4)*f2*h**3*x0*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi + (3/4)*f2*h**2*x0*x2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi - 0.75*f2*h*x0*(h + x2)**2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi - 1/8*f2*h*x0*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi + (1/8)*f2*h*x0*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + (1/8)*f2*x0*x2*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi - 1/8*f2*x0*x2*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi)

lambda x0,x1,x2,f0,f1,f2,h: (0.75*f0*h*x0*x1*x2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi + (1/8)*f0*x0*x1*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi - 1/8*f0*x0*x1*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + 0.75*f1*h*x1**2*x2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi - 1/4*f1*h*x2*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + (1/8)*f1*x1**2*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi - 1/8*f1*x1**2*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + (1/8)*f1*(x0**2 + x1**2 + (h - x2)**2)**(-0.5)/math.pi - 1/8*f1*(x0**2 + x1**2 + (h + x2)**2)**(-0.5)/math.pi + (3/4)*f2*h**3*x1*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi + (3/4)*f2*h**2*x1*x2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi - 0.75*f2*h*x1*(h + x2)**2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi - 1/8*f2*h*x1*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi + (1/8)*f2*h*x1*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + (1/8)*f2*x1*x2*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi - 1/8*f2*x1*x2*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi)

lambda x0,x1,x2,f0,f1,f2,h: (-3/4*f0*h**3*x0*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi - 3/4*f0*h**2*x0*x2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi + 0.75*f0*h*x0*(h + x2)**2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi - 1/8*f0*h*x0*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi + 0.125*f0*h*x0*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + (1/8)*f0*x0*x2*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi - 1/8*f0*x0*x2*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi - 3/4*f1*h**3*x1*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi - 3/4*f1*h**2*x1*x2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi + 0.75*f1*h*x1*(h + x2)**2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi - 1/8*f1*h*x1*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi + 0.125*f1*h*x1*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + (1/8)*f1*x1*x2*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi - 1/8*f1*x1*x2*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi - 0.75*f2*h*x2*(h + x2)**2*(x0**2 + x1**2 + (h + x2)**2)**(-2.5)/math.pi + 0.25*f2*h*x2*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + (1/8)*f2*(h - x2)**2*(x0**2 + x1**2 + (h - x2)**2)**(-1.5)/math.pi - 1/8*f2*(h + x2)**2*(x0**2 + x1**2 + (h + x2)**2)**(-1.5)/math.pi + (1/8)*f2*(x0**2 + x1**2 + (h - x2)**2)**(-0.5)/math.pi - 1/8*f2*(x0**2 + x1**2 + (h + x2)**2)**(-0.5)/math.pi)



In [26]:

    
t1 = []
lam_u0 = lambdify((x0, x1, x2, f0, f1, f2, h), u0)
for _ in range(1000):
    t1.append(lam_u0(*np.random.sample(2), 0, 1, 1, 1, 1))
print(np.sum(t1))

0.0



In [16]:

    
u0, u1, u2 = simplify(Mij * F)
display(u0)
display(u1)
display(u2)









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-16-b35418f0c1bc> in <module>
----> 1 u0, u1, u2 = simplify(Mij * F)
      2 display(u0)
      3 display(u1)
      4 display(u2)

NameError: name 'Mij' is not defined



In [14]:

    
type(divU)









    Out[14]:





sympy.core.add.Add