[Python] Notebook-12 Notes about SymPy

這一篇仍然不是典型的上課筆記，
是我自己研究SymPy時的心得紀錄。
從11/12/2017開始寫，
到現在已經一個多月了。
這段時間除了看上面這個SymPy的網站，
也複習電子學、微分方程和拉普拉斯轉換(Laplace Transform)。

說到微分方程就必須要推薦材料工程學系的徐雍鎣老師教授的工程數學一，

清楚又有結構性，
很多人都說教得好。

我自己上課時把筆記再整理過一次並加上書籤，
有上傳到部落格的分享資料夾裡，

就在部落格首頁的右手邊而已；

也有分享在Youtube的第一週課程影片影片下方，
歡迎自行下載。

好，
回到這篇文章的正題—SymPy。
這一篇文章主要內容來自官方索引（？！就是本文第一個超連結啦！），
只是加上一點例子讓大家看起來比較方便。
內容包含開根號的符號表示法和近似值、多項式的展開和取公因式、不定積分、定積分和微分、指數運算、迭代運算(lambdify)、微分方程求解與繪圖。
以下附上運算與顯示結果給大家參考。

看看下面的微分方程例題二，
就會發現一個特殊的函數LambertW，
後面的例題還會出現大O符號。
這樣的解其實不太可能用手算得到。
因為這樣的微分方程通常沒有解析解（公式解），
所以電腦幫忙求出來的多是近似解，
這一點和WolframAlpha就大不相同。
至於WolframAlpha的使用，
我就不多談，
菇狗是你的好幫手。

大家可以看一下下面的微分方程範例第八題，
因為求出來的解是複數解，
沒有辦法在二維的平面作圖，
所以下面繪圖的部分我就跳過去了。

以下就是程式碼的部分。

from __future__ import division #It must be at the begining.
from sympy import * #This imports all the functions and classes from SymPy.
x, y, z, t, C1 = symbols('x y z t C1')
k, m, n = symbols('k m n', integer = True)
f, g, h, O = symbols('f g h O', cls = Function)
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
The above five rows are automatically executed once you clicked the check 
box of "Use symbolic math." But if you do not check the box, remember to type
these commands to make it work in the way you expect. More information could 
be found here: http://blog.csdn.net/u012675539/article/details/46981305
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

# -*- coding: utf-8 -*-
"""
Created on Fri Dec  8 11:36:46 2017
@author: ShihHsing Chen

The following note are cut from "SymPy Tutorial." You can go to the webpage 
for more information. http://docs.sympy.org/latest/tutorial/index.html
"""

#from sympy import symbols
#from sympy import expand, factor
#from sympy import Integral, preview
#from sympy import init_printing
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
The function of the four rows above is replaced by the 2nd row of this file:
from sympy import *
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
import matplotlib.pyplot as plt
import time #I just want to know how much time does a run would take.
start = time.time()

#"Color" should be titlized, and backcolor is also available. This function 
#is quite useful as you work under a dark theme. It would give you a clear 
#view of your "display()."
init_printing(use_latex = True, forecolor = "Green") 

#Symbolic computation with sympy
import sympy
sqrt8 = sqrt(8)
print("symbolic expression of square root of 8:")
display(sqrt8)

#To evaluate a numerical expression into a floating point number, use 
#.evalf(n). n is the precision you want, and the default value is 15.
num_sqrt8 = sqrt8.evalf()
print("numerical expression of square root of 8:")
display(num_sqrt8)

#In sympy, variables are defined before using. (Just like other package.) 
expr = x**2 + 2*x*y + y**2 + 100
print("original polynomial:")
display(expr)

#In sympy, there are functions to use for expanded form or factored form.
factored_expr = factor(x*expr)
expanded_expr = expand(x*expr)
print("original polynomial times x in factored form:")
display(factored_expr)
print("original polynomial times x in expanded form:")
display(expanded_expr)

#SymPy could also help you solving Calculus. Note that SymPy does not include 
#the constant of integration. You could read more here:
"""https://goo.gl/QrgEVU"""
print("If you have problem solving an indefinite integration like this one,")
a = Integral(cos(x)*exp(x), x)
display(a)
print("the method .doit will execute the calculation for you.")
display(Eq(a, a.doit())) #It will not pop out without display.
                         #Eq(a, b) is used to create symbolic equalities. 
                         #Read more below.
print()

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
                               Equals signs

In SymPy, == represents EXACT structural equality testing. This means that 
a == b means that we are asking if a=b. We always get a bool as the result 
of ==. There is a separate object, called Eq, which can be used to create 
symbolic equalities.

There is also a method called equals that tests if two expressions are equal 
by evaluating them numerically at random points. If a equals b, the command
a.equals(b) will give you a "True."
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
                                  ^ and /

You may have noticed that we have been using ** for exponentiation instead 
of the standard ^. That’s because SymPy follows Python’s conventions. In 
Python, ^ represents logical exclusive or. SymPy follows this convention.

Whenever you combine a SymPy object and a SymPy object, or a SymPy object 
and a Python object, you get a SymPy object, but whenever you combine two 
Python objects, SymPy never comes into play, and so you get a Python object.

This is usually not a big deal. Python ints work much the same as SymPy 
Integers, but there is one important exception: division. In SymPy, the 
division of two Integers gives a Rational. But in Python / represents either 
integer division or floating point division, depending on whether you are in 
Python 2 or Python 3, and depending on whether or not you have run 
"from __future__ import division."

To avoid this, we can construct the rational object explicitly Rational(a, b).
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
                               Substitution

One of the most common things you might want to do with a mathematical 
expression is substitution. Substitution replaces all instances of something
in an expression with something else. It is done using the .subs() method. 
Substitution is usually done for one of two reasons:
 
1. Evaluating an expression at a point. For example, if our expression is 
cos(x) + 1 and we want to evaluate it at the point x = 0, so that we get 
cos(0) + 1, which is 2.

2. Replacing a subexpression with another subexpression. There are two 
reasons we might want to do this. The first is if we are trying to build an 
expression that has some symmetry, such as x**x**x**x. To build this, we 
might start with x**y, and replace y with x**y. We would then get x**(x**y). 
If we replaced y in this new expression with x**x, we would get 
x**(x**(x**x)), the desired expression.

The second is if we want to perform a very controlled simplification, or 
perhaps a simplification that SymPy is otherwise unable to do. For example, 
say we have sin(2x)+cos(2x), and we want to replace sin⁡(2x) with 
2sin⁡(x)cos⁡(x). As we will learn later, the function expand_trig does this. 
However, this function will also expand cos⁡(2x), which we may not want. 
While there are ways to perform such precise simplification, and we will 
learn some of them in the advanced expression manipulation section, an easy 
way is to just replace sin(2x) with 2sin(x)cos(x).

There are two important things to note about subs. First, it returns a new 
expression. SymPy objects are immutable. That means that subs does not 
modify it in-place.

To perform multiple substitutions at once, pass a list of 
([(a, b), (c, d), (e, f)]) pairs to subs. It is often useful to combine this 
with a list comprehension to do a large set of similar replacements all at 
once.

To numerically evaluate an expression with a Symbol at a point, we might 
use subs followed by evalf, but it is more efficient and numerically stable 
to pass the substitution to evalf using .evalf(subs={x: a}).

Sometimes there are roundoff errors smaller than the desired precision that 
remain after an expression is evaluated. Such numbers can be removed at the 
user’s discretion by setting the chop flag to True. Like .evalf(chop = True).
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

print("x to the power of y:")
expr = x**y
display(expr)
print("substitution of y with x to the power of y:")
expr = expr.subs(y, x**y)
display(expr)
print("substitution of y with x to the power of x:")
expr = expr.subs(y, x**x)
display(expr)

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
                   Evaluation at Multiple Points: Lambdify

subs and evalf are good if you want to do simple evaluation, but if you 
intend to evaluate an expression at many points, there are more efficient 
ways. For example, if you wanted to evaluate an expression at a thousand 
points, using SymPy would be far slower than it needs to be, especially if 
you only care about machine precision. Instead, you should use libraries 
like NumPy.

The easiest way to convert a SymPy expression to an expression that can be 
numerically evaluated is to use the lambdify function. lambdify acts like a 
lambda function, except it converts the SymPy names to the names of the 
given numerical library, usually NumPy. We take a look at the following 
example.
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
import numpy as np
print()
print("Let's evaluate the value of sinx with x from 0 to 9.")
array = np.arange(10) 
print("the array of evaluation points:", "\n", array)
expr_to_convert = sin(x)
expr_converted1 = lambdify(x, expr_to_convert, "numpy")
print("the array of evaluation results:", "\n", expr_converted1(array))

#The topic followed is "Simplify," and I am going to leave it for you to 
#explore. The content covers (1) Polynomial/Rational function simplification 
#(2) Trigonometric simplification (3) Powers (4) Exponentials and logarithms 
#(5) Special functions.
"""https://goo.gl/vhrYTL"""

#The topic followed is "Calculus," and I am going to leave it for you. The 
#content covers (1) Derivatives (2) Integrals (3) Limits (4) Series expansions
#(5) Finite differences. You could discover the topic with your
#curiousity and passion.
"""https://goo.gl/QrgEVU"""

#The topic followed is "Solvers," and I am going to skip it except the part of
#differential equations below. To solve differential equations, we use dsolve.
#dsolve(eq, func=None, hint='default', simplify=True, ics=None, xi=None, 
#eta=None, x0=0, n=6, **kwargs)
"""https://goo.gl/9SZ1Fh"""
#Modules Reference of ODE
"""https://goo.gl/5b3VxH"""
#Modules Reference of PDE
"""https://goo.gl/njHEuo"""

#I have done problem sets in the first chapter of Prof. Yung-Jung Hsu's 
#teaching materials. Check it out!
"""https://goo.gl/1Dsn7r"""
#Besides, I am inspired to know how to plot the differential equations by 
#Bill Bell's answer here.
"""https://goo.gl/z1j8pk"""
print("\n")
print("Differential equation example one: Separable Equation")
diffeq1 = Eq(f(x).diff(x), f(x)/(1+x))
display(diffeq1)
print("Solution:")
gen_sol1 = dsolve(diffeq1, f(x))
display(gen_sol1)
plot1 = gen_sol1.rhs.subs(C1, 1)
plot(plot1, (x, -100, 100), title = "Solution with C1 = 1", autoscale= True)

print("Differential equation example two: Separable Equation")
diffeq2 = Eq(f(x).diff(x), (x*f(x)+3*x-f(x)-3)/(x*f(x)-2*x+4*f(x)-8))
display(diffeq2)
print("Solution:")
gen_sol2 = dsolve(diffeq2, f(x))
display(gen_sol2)
plot2 = gen_sol2.rhs.subs(C1, 1)
plot(plot2, (x, -100, 100), title = "Solution with C1 = 1", autoscale= True)

print("Differential equation example three: Separable Equation")
diffeq3 = Eq(f(x).diff(x), f(x)/x+1)
display(diffeq3)
print("Solution with C1 = 1:")
gen_sol3 = dsolve(diffeq3, f(x))
display(gen_sol3)
plot3 = gen_sol3.rhs.subs(C1, 1)
plot(plot3, (x, -100, 100), title = "Solution with C1 = 1", autoscale= True)

print("Differential equation example four: Separable Equation")
diffeq4 = Eq(f(x).diff(x), (x-f(x))/(x+f(x)))
display(diffeq4)
print("Solution:")
gen_sol4 = dsolve(diffeq4, f(x))
display(gen_sol4)
plot4a = gen_sol4[0].rhs.subs(C1, 1)
plot4b = gen_sol4[1].rhs.subs(C1, 1)
plot(plot4a, (x, -100, 100), title = "Solution A with C1 = 1", autoscale= True)
plot(plot4b, (x, -100, 100), title = "Solution B with C1 = 1", autoscale= True)

print("Differential equation example five: Separable Equation")
diffeq5 = Eq(f(x).diff(x), tan(x+f(x))*tan(x+f(x)))
display(diffeq5)
print("Solution:")
gen_sol5 = dsolve(diffeq5, f(x))
display(gen_sol5)
plot5 = gen_sol5.rhs.subs(C1, 1)
plot5 = plot5 - O(x**6)
plot(plot5, (x, -100, 100), title = "Solution with C1 = 1", autoscale= True)
#Example five would generate an approximate solution, which is unlikely 
#done by hand-calculations. As a result, I suggest readers use the 
#WolframAlpha, which would give you a equation that you could confirm with
#hand-calculations.
"""https://goo.gl/T4JYP8"""

print("Differential equation example six: Exact equation")
diffeq6 = Eq(f(x).diff(x), (sin(f(x))-f(x)*sin(x))/(f(x)-cos(x)-x*cos(f(x))))
display(diffeq6)
print("Solution:")
gen_sol6 = dsolve(diffeq6, f(x))
display(gen_sol6)
plot6 = gen_sol6.rhs.subs(C1, 1)
plot6 = plot6 - O(x**6)
plot(plot6, (x, -100, 100), title = "Solution with C1 = 1", autoscale= True)

print("Differential equation example seven: Modified exact equation")
diffeq7 = Eq(f(x).diff(x), (-6*x*f(x)/(4*f(x)-9*x**2)))
display(diffeq7)
print("Solution:")
gen_sol7 = dsolve(diffeq7, f(x))
display(gen_sol7)
plot7 = gen_sol7.rhs.subs(C1, 1)
plot7 = plot7 - O(x**6)
plot(plot7, (x, -100, 100), title = "Solution with C1 = 1", autoscale= True)

print("Differential equation example eight: Linear Differential Equation")
diffeq8 = Eq(f(x).diff(x), -f(x)*tan(x)-sin(2*x))
display(diffeq8)
print("Solution:")
gen_sol8 = dsolve(diffeq8, f(x))
display(gen_sol8)

print("Differential equation example nine: Bernoulli’s Equation")
diffeq9 = Eq(f(x).diff(x), f(x)**3-2*f(x)/x)
display(diffeq9)
print("Solution:")
gen_sol9 = dsolve(diffeq9, f(x))
display(gen_sol9)
plot9 = gen_sol9.rhs.subs(C1, 1)
plot(plot9, (x, -100, 100), title = "Solution with C1 = 1", autoscale= True)

#The topic followed is "Matrices," and I am going to skip it. If you need to 
#use matrices at higher dimensions, search numpy.ndarray at my blog.
"""https://goo.gl/9SZ1Fh"""

end = time.time()
elapsed = end - start
print("Time taken: ", elapsed, "seconds.")