Sympy：求解球面与平面的交线，并画出图

• 用Sympy来解空间解析几何
• 先上全部代码，后面有详细说明的Jupyter notebook版本

# -*- coding: utf-8 -*-
"""
Created on Thu Aug 30 20:11:37 2018

@author: phart
"""
import sympy as sp
import numpy as np
from sympy import *
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as p3d
from sympy.abc import alpha, beta, gamma, delta, mu, sigma, epsilon
from IPython.core.interactiveshell import InteractiveShell

sp.init_printing(use_unicode=True)
InteractiveShell.ast_node_interactivity = 'all'

x, y, z = symbols('x, y, z')
i, n = symbols('i, n')

eq = [x**2 + y**2 + z**2 - 1,
      x - z]
s1, s2 = nonlinsolve(eq, (x, y))
x1, y1 = s1
x2, y2 = s2

zm = np.sqrt(1/2) - 0.0005         # 由于上面的解是z的函数，所以我们以z取样，
zl = np.linspace(-zm, zm, 50)      # 这里减0.0005是由于精度问题，不减的话会有虚数出现
xil = []                           # 因为上面方程的解是一个空间圆线，所以x, y都分成两部分
xir = []
yir = []
yil = []
for zi in zl:                      # 带入z, 求解相应的x, y
    xil.append(x1.subs({z:zi}).doit().evalf(4))
    xir.append(x2.subs({z:zi}).doit().evalf(4))
    yil.append(y1.subs({z:zi}).doit().evalf(4))
    yir.append(y2.subs({z:zi}).doit().evalf(4))

fig = plt.figure()
ax = p3d.Axes3D(fig)

ax.plot(xil, yil, zl, color='g')
ax.plot(xir, yir, zl, color='g')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

ax.view_init(30, 120)

plt.show()

• 效果如图

• 下面是notebook版本的，还是粘贴图片比较方便 ~~

  
  
  
  
• 最后一行ax.view_init(30, 120)是用来旋转坐标轴，达到比较好的输出视角，第一个参数30代表俯视的角度，第二个参数是代表以z轴为轴旋转的角度。
• 至于最后一个图中等号前多出的下划线是为了输出好看一些，因为在第一个图中设置了所有单行变量都会输出，所以就把作图返回的对象赋值给下划线，可以在notebook中把下划线和等号删除试验一下。