блог JOHN_16: Python: модуль SymPy. Пример символьных вычислений.

SymPy - это обширный Python модуль реализующий символные вычисления. В качестве примера будет разобрана реальная математическая задача.

Суть задачи:

Дана окружность радиусом r с координатами центра в точке (x₀, y₀). Также дана точка с координатами (x₁, y₁). Необходимо провести лучи из центра окружности, проходящие через точку пересечения окружности и касательных. Все элементы задачи находятся на одной плоскости.

Задача составная, поэтому разделим ее на логически отделенные части, решающие наиболее важные асппекты: 1) нахождение касательных 2)нахождение точек пересечения касательной и окружности. После каждой части будет приведен фрагмент кода, подкрепленный комментариями.

Часть 1. Нахождение касательных

Для начала решим задачу в общем виде. Запишем математическое выражение окружности и касательной в виде системы уравнений:

На два уравнения приходится три неизвестных (x,y,k), с математической точки зрения система имеет множество решений, однако далее мы введем условие, которое разрешит данное положение. Из уравнения прямой выразим y:

И подставим в уравнение окружности, далее, обобщив, получим квадратное уравнение:

По определению касательная и окружность имеют только одну точку пересечения. Это возможно только в случае, если квадратное уравнение имеет 1 корень. По определению квадратного уравнения оно имеет 1 корень в случае, если дискриминант равен нулю:

Решая полученное уравнение получаем несколько значений k₁, k₂. Подставляем значения в общее уравнение прямой и получаем пару уравнений касательных проведенных из точки (x₁, y₁) к окружности

Для демонстрации численного решения воспользуемся следующими начальными данными: окружность радиусом 4 с центром в точке (-5, 1) и точка касательной с координатами (10, 15).

Код с комментариями:

# -*- coding: utf-8 -*-

from sympy import *
import numpy as np

# точка центра окружности
x0 = -5
y0 = 1

# радиус окружности
r = 4

# точка вершины прямой
x1 = 10
y1 = 15

# Определяем переменные
x, y, k = symbols('x y k')

# Выражаем y из уравнения касательной
exp_2 = k * (x - x1) - (y - y1)

# решаем полученное выражение относительно y
# функция solve всегда возвращает список решений, но поскольку мы заранее знаем, что 
# решений может быть только одно, явно извлекаем нулевой элемент
ny = solve(exp_2, y)[0]
print 'y=', ny

# Выражаем x из уравнений окружности, подставляя ранее полученное
# выражение y
exp = (x - x0) ** 2 + (ny - y0) ** 2 - r ** 2

# приводим выражение к общему виду относительно переменной x функцией collect,
# метод cancel позволяет "упростить" выражение, преобразовав подобные элементы
ne = collect(exp.cancel(), x)

# преобразуем выражения в объект полинома
po = poly_from_expr(ne, x)[0]

# выделяем коэфициенты квадратного уравнения
alpha = po.coeff_monomial(x**2)
beta = po.coeff_monomial(x)
gamma = po.coeff_monomial(1)
print 'alpha =', alpha
print 'beta =', beta
print 'gamma =', gamma

# Касательная и окружность пересекаются в одной точке, поэтому уравнение
# должно иметь только 1 корень, это происходит в случае когда дискриминант
# равен нулю. Решаем уравнение и находим значения углового коэффициента
D_exp = beta ** 2 - 4 * alpha * gamma
res = solve(D_exp, k)

# evalf метод вычисляет выражение в численной форме
# здесь мы его используем что бы получить более наглядное представление 
print 'k =', res, '~~>', [_.evalf() for _ in res]

# Подставляем полученные значения k и получаем уравнения касательных.
tangents_lines = []
for item in res:
    # метод subs осуществляет подстановку значений в переменные и возвращает 
    # новое выражение
    tangents_lines.append(ny.subs({k:item}))
    print 'y_tangent =', tangents_lines[-1] , '~~>', tangents_lines[-1].evalf()

Вывод данного фрагмента кода будет следующим:

y= k*x - 10*k + 15
alpha = k**2 + 1
beta = -20*k**2 + 28*k + 10
gamma = 100*k**2 - 280*k + 205
k = [-36*sqrt(5)/209 + 210/209, 36*sqrt(5)/209 + 210/209] ~~> [0.619624654593338, 1.38994472339709]
y_tangent = x*(-36*sqrt(5)/209 + 210/209) + 360*sqrt(5)/209 + 1035/209 ~~> 0.619624654593337*x + 8.80375345406662
y_tangent = x*(36*sqrt(5)/209 + 210/209) - 360*sqrt(5)/209 + 1035/209 ~~> 1.38994472339709*x + 1.10055276602907

Часть 2. Поиск точек пересечений касательных и окружности

Библиотека SymPy для задач аналитической геометрии предоставляет геометрические объекты и различного рода функции для работы с ними. В продолжение прошлой части кода, зададим такие объекты и найдем их точки пересечений, которыми воспользуемся для определения других геометрических объектов.

# список в котором будут храниться точки пересечений
intersections = []
# задаем объект окружности
circle = Circle(Point(x0, y0), r)
# задаем точку из которой будем строить касательную
main_point = Point(x1, y1)

# для каждой касательной ищем пересечение с окружностью
for tangent in tangents_lines:
    
    # задаем вторую точку касательной
    new_point = Point(x0, tangent.subs({x:x0}))
    
    # если вдруг так получилось что точки совпали, то зададим небольшое смещение
    # и пересчитаем координаты точки
    if main_point == new_point:
        new_point = Point(x0 + 1, tangent.subs({x:x0 + 1}))
    
    # задаем объект линию, проходящую через заданные ранее точки
    line = Line(main_point, new_point)
    
    # находим пересечение двух геометрических объектов с помощью 
    # функции intersection, которая вернет список точек пересечений
    tmp = intersection(circle, line)
    
    # если список tmp не пуст, значит пересечения были найдены
    if tmp:
        # если длина списка равна единице, то значит все хорошо
        # для данного случая - мы нашли одну единственну точку пересечения
        # что соответсвует определению касательной к окружности
        if len(tmp) == 1:
            intersections.append(tmp[0])
        
        # если вдруг нашлось более одного пересечения, то в каком то уголке вселенной 
        # нерадивый ученик разделил на ноль, извлек корень из минус единицы 
        # и определил последнее число Пи после запятой. 
        else:
     print 'I don\'t know how it was, but this is your points: {}'.format(tmp)
    else:
        print 'Circle "{}" and line "{}" have no intersection'.format(circle, line)

#
#
# Строим лучи с вершиной в центре окружности, проходящие через точки пересечения с касательной
lines = []
for inter in intersections:
    lines.append(Ray(Point(x0,y0), inter))
    print 'Ray', lines[-1], 'pretty views points of ray', [_.evalf() for _ in lines[-1].args]

Вывод данного фрагмента кода будет следующим:

Ray Ray(Point(-5, 1), Point(-1865/421 - 504*sqrt(5)/421, 645/421 + 540*sqrt(5)/421)) pretty views points of ray [Point(-5.00000000000000, 1.00000000000000), Point(-7.10683672365771, 4.40018220391897)]
Ray Ray(Point(-5, 1), Point(-1865/421 + 504*sqrt(5)/421, -540*sqrt(5)/421 + 645/421)) pretty views points of ray [Point(-5.00000000000000, 1.00000000000000), Point(-1.75302075852757, -1.33604918729189)]

Заключение

На рисунке ниже графически представлены некоторые результаты вычислений. Черным цветом нарисованы исходные положения - окружность и точка. Красными оттенками нарисованы две касательные, построенные по вычиленным уравнениям. Зеленым цветом отмечены лучи, исходящие из центра окружности, и, проходящие через точки пересечений касательной и окружности, которые были расчитаны во второй части этого сообщения.

Поставленная задача успешно решена с помощью возможностей предоставляемой библиотекой SymPy.