I'm addicted to formulas in python

What I'm learning now is so interesting when I try various things with python, and I forgot to write an article. You don't have to write python in code for substitution, factorization, differentiation, etc. The library is too convenient. I am currently studying elementary chaos theory, multiple regression analysis of multivariate analysis, and multivariate normal distribution as elective courses at the university. I will visualize it because I can not understand anything even if I look at the symbols that I do not understand well in the classroom. I don't understand the theory at all, but if you look at it, it may be a chance to gain some understanding. Chaos theory using multivariate analysis and differential equations is too difficult, so I'm stumbling from the beginning, so this time I searched variously on the net from the beginning to the analysis, which is ideal I've been searching for the code. To say the difficult words, I don't understand anything ... I didn't write this article myself .... Difficult ..... Someday I want to understand the theory of multivariate analysis and code it. No, why did you get so motivated when you didn't understand anything?

Program code basic library

I imported the function with sympy in python. Bruise the library. It's so convenient that the derivative is displayed so intuitively ... As you can see from the code, the function that uses differentiation is diff ().

import sympy as sym
from sympy import *      #Imoprt all features from the sympy library
sym.init_printing()      #Format output

x=Symbol('x')                  #letter'x'Is defined as the variable x
y=Symbol('y')                 #letter'y'Is defined as the variable y


#differential
a=diff(sin(x),x)  #sin(x)Differentiation of
b=diff(exp(x),x)   # e^Differentiation of x

#Higher derivative
c=diff(x**4+x**3,x,2)  # x^4+x^2nd derivative of 3


#Partial differential
d=diff(x**2+x*y+2*y**2,x)  # x^2+xy+2y^Partial derivative of 2 by x
e=sym.integrate(x+2*y,(x,0,2),(y,2,3))
print('a: {0}'.format(a))
print('b: {0}'.format(b))
print('c: {0}'.format(c))
print('d: {0}'.format(d))
print('e: {0}'.format(e))

From top to bottom

cos(x)\\
e^x\\
6x(2x + 1)\\
2x+y\\
24

Graph display of two-variable function


{3x^2+6xy+3y^2+x^3+y^3}

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
def func1(x, y):
    return 3*x**2+6*x*y+3*y**2+x**3+y**3
x = np.arange(-10.0, 5.0, 0.1)
y = np.arange(-10.0, 5.0, 0.1)
X, Y = np.meshgrid(x, y)
Z = func1(X, Y)
fig = plt.figure()
ax = Axes3D(fig)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("f(x, y)")
ax.plot_wireframe(X, Y, Z)
plt.show()

It seems that the point of the maximum value is at the place of (-4, -4) when the troublesome stop point and the maximum point are calculated by hand. It's hard to understand. I wonder if it matches

Find the local maximum value

from sympy import *


x = Symbol('x')

y = Symbol('y')



f =3*x**2+6*x*y+3*y**2+x**3+y**3




f_x = Derivative(f,x).doit() #First-order partial differential with x

f_y = Derivative(f,y).doit() #First-order partial differential with y



stationary_points = solve([f_x,f_y],[x,y]) #Function stop point

stationary_points = list(map(lambda x:simplify(x),stationary_points)) #Simplification



f_xx = Derivative(f, x, 2).doit() #Second-order partial derivative with x

f_yy = Derivative(f, y, 2).doit() #Second-order partial derivative with y

f_xy = Derivative(f_x,y).doit() #Partial differentiation with x and then partial differentiation with y

f_yx = Derivative(f_y,x).doit() #Partial differentiation with y and then partial differentiation with x



#Hessian matrix

hesse = Matrix([

[f_xx,f_xy],

[f_yx,f_yy]

])



#Determinant of Hessian matrix

det_hessian = hesse.det()



#Output for the time being

print('f_x: {0}'.format(f_x))

print('f_y: {0}'.format(f_y))

print('f_xx: {0}'.format(f_xx))

print('f_yy: {0}'.format(f_yy))

print('f_xy: {0}'.format(f_xy))

print('f_yx: {0}'.format(f_yx))

print('All stops{0}'.format(stationary_points))



#Judgment

for i in range(len(stationary_points)):

    sign_det_hessian = det_hessian.subs([(x,stationary_points[i][0]),(y,stationary_points[i][1])]).evalf()

    if sign_det_hessian > 0:

        sign_f_xx = f_xx.subs([(x,stationary_points[i][0]),(y,stationary_points[i][1])]).evalf()

        if sign_f_xx > 0:

            print('Stop point:{0}Is f(x,y)Minimal point'.format([stationary_points[i][0],stationary_points[i][1]]))

        elif sign_f_xx < 0:

            print('Stop point:{0}Is f(x,y)Maxima'.format([stationary_points[i][0],stationary_points[i][1]]))

        else:

            print('Cannot be determined by this method')

    elif sign_det_hessian < 0:

        print('Stop point:{0}Is f(x,y)Saddle point'.format([stationary_points[i][0],stationary_points[i][1]]))

There was, but I don't understand the meaning. You can't even see it with your eyes. Only the image of a handkerchief flying in the air comes to mind ...

Make sure you are at the limit

I checked to see if it was the limit. Using polar coordinates, we used x = rcosθ y = 2rsinθ Change the polar coordinate value according to the denominator value and the situation, create a trigonometric function formula, etc., and try about a fraction. By the way, the convergence value is 2.


\lim_{(x,y)\to (0,0)}\frac{x^3-8y^3+2x^2+8y^2}{x^2+4y^2}

import sympy as sym
import math
sym.init_printing()
Pi = sym.S.Pi #Pi
E = sym.S.Exp1 #The bottom of the natural logarithm
I = sym.S.ImaginaryUnit #Imaginary unit
oo = sym.oo #Infinity

#Definition of variables to use(All lowercase letters are symbols)
(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z) = sym.symbols('a b c d e f g h i j k l m n o p q r s t u v w x y z')

import numpy as np
import matplotlib.pyplot as plt
c=(x**3-8*y**3+2*x**2+8*y**2)/(x**2+4*y**2)
d=sym.limit(c,x,0)
print('{0}'.format(d))
dd=sym.limit(d,y,0)
e=sym.limit(c,y,0)
ee=sym.limit(e,x,0)
print('{0}'.format(e))
if dd == ee :
    print('The limit of a function exists')
    f=c.subs([(x,X),(y,1/2*Y)]).trigsimp() #Need simplification(You should do it line by line)
    g=np.abs(sym.limit(f-ee,r,0))
    print('{0}'.format(g))
else :
    print('There is no limit of function')

This execution result shows that the limit of the function exists safely, but when the limit does not exist, we have not verified how it is executed. I will do my best.

Display of execution result of partial differentiation


(x^2+xy+y^2)\log(x^2+y^2)

I tried to partially differentiate x and y respectively.

from sympy import *
import sympy as sym
import math
from IPython.display import display, Latex
x=Symbol('x')                  #letter'x'Is defined as the variable x
y=Symbol('y')                 #letter'y'Is defined as the variable y

f = (x ** 2 + x * y + y ** 2) * sym.log(x ** 2 + y ** 2)

display(f)

dx = sym.diff(f, x).subs({y: 0})
dy = sym.diff(f, y).subs({x: 0})

display(dx)
sym.plot(dx, (x, -5, 5))

display(dy)
sym.plot(dy, (y, -5, 5))

When I accepted the question at teratail, I was able to get the answer. Thank you very much. If x and y are partially differentiated, both pass through the origin, so partial differentiation is possible at the origin. Easy to understand.

Impressions

It was really convenient for python. Compared to other languages, it has a lot of libraries and is really amazing. Next, I would like to add a theory. This time too, it was just an introduction. I want to do something more applied. Do you use this code for troublesome calculations and checkups?

Reference sites and questions

https://qiita.com/tibigame/items/aebbac176d9bbdaf3d15 https://qiita.com/hiroyuki_mrp/items/d373b951e216f62c4957 https://qiita.com/PlanetMeron/items/63ac58898541cbe81ada https://short4010.hatenablog.com/entry/2019/05/12/142731 It's my account https://teratail.com/questions/283867

I tried Jacobian and partial differential with python