Введение в научный Python-1. Введение в научный Python
 Скачать 6.2 Mb.
|
|
Пример. Лента Мебиуса – простейшая неориентируемая поверхность. Она получается движением и вращением отрезка прямой вдоль замкнутой пространственной кривой. Пусть эта кривая будет окружностью радиуса R, а n обозначает количество полуоборотов отрезка при обходе кривой. Вдоль направляющей окружности движется центр отрезка. Уравнение поверхности можно записать в параметрическом виде u u n v R v u x cos 2 cos , u u n v R v u y sin 2 cos , 2 sin , u n v v u z Пусть H обозначает ширину ленты (длину отрезка). При нечетном n получается неориентированная поверхность (лента Мебиуса). 169 Следующий код строит анимацию движения отрезка и заметаемую им поверхность. import numpy as np import matplotlib.pyplot as plt from mpl_too lkits.mplot3d import Axes3D import matplotlib.animation as animation # координатные функции xs=lambda u,v: (R+v*np.cos(n*u/2))*np.cos(u) ys=lambda u,v: (R+v*np.cos(n*u/2))*np.sin(u) zs=lambda u,v: v*np.sin(n*u/2) def initpict(): global u u=np.linspace(0,du,1) # очистка вектора u перед стартовым кадром def redraw(i): # рисование i-го кадра global u u=np.append(u,(i+1)*du) # добавление точки к вектору u U,V=np.meshgrid(u,v) X=xs(U,V) Y=ys(U,V) Z=zs(U,V) ax.clear() # после очистки восстанавливаем видимые размеры графической области ax.set_xlim3d(xyLeft, xyRight) ax.set_ylim3d(xyLeft, xyRight) ax.set_zlim3d(-H,H) ax.plot(xc,yc,zc,color='brown',linewidth=3) # направляющая окружность ax.plot_surface(X, Y, Z, rstride=1,cstride=1,color='0.5',alpha=0.4) R=4;H=2;n=1; NumFrames=40 # количество кадров du=2*np.pi/NumFrames xyRight=R+1; xyLeft= -xyRight; fig = plt.figure(facecolor='white') fig.set_dpi(100) ax=Axes3D(fig,xlim=(xyLeft,xyRight),ylim=(xyLeft,xyRight),zlim=( -H,H)) # координаты точек направляющей окружности одинаковые для всех кадров t=np.linspace(0,2*np.pi,100) vt=np.zeros(np.size(t)) xc=xs(t,vt) yc=ys(t,vt) zc=zs(t,vt) # массив параметра v поверхности одинаков для всех моментов времени v=np.linspace( -H/2,H/2,11) anim =animation.FuncAnimation(fig, redraw, init_func=initpict, 170 frames=NumFrames, interval=20, repeat=True,blit=False) plt.show() На следующем рисунке показано несколько кадров анимации в последовательные моменты времени (при значении параметра n=1). Особенность этого кода состоит в том, что в каждом кадре, создаваемом функцией redraw, добавляется только одна «вертикальная полоска» поверхности. Для этого перед началом анимации функция инициализации initpict() очищает вектор u, который в результате будет содержать только нулевое значение. Функция redraw в каждом новом кадре к вектору u добавляет один элемент(i+1)·du, больший предыдущего на величину NumFrames / 2 u d . Следующие четыре команды создают массивы U,V значений параметров и массивы X, Y, Z координат точек участка поверхности ленты. Но прежде, чем ленту рисовать, инструкция ax.clear() очищает графическую область, а последующие три инструкции восстанавливают для нового изображения прежние размеры видимой области. После этого рисуется направляющая окружность, а затем рисуется участок ленты u d 1 u 0 i Если очистку графического окна не делать, то в графической области после создания i-го кадра будет нарисовано i налагающихся друг на друга кусков ленты разной длины. Сверху будет нарисован самый длинный, но в памяти будут присутствовать все. Это значительно «утяжеляет» рисунок, который будет сложнее «крутить». Первый кадр анимации строится инструкцией redraw(0), которая вызывается из функции FuncAnimation. Она рисует одну полоску u d u 0 Последующий вызов redraw(1) рисует две полоски ( u d 2 u 0 ) и т.д. Протестируйте работу программы при n=2 и n=3. ■ Иногда кадры анимации содержат изображения матриц, взятых из файла или сгенерированных в коде программы. Ниже приведен код, создающий анимацию изображений матриц, которые «рисуются» функцией imshow. Пример. import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation fig = plt.figure(facecolor='white') f=lambda x,y: np.sin(2*x)+np.sin(2*y) x = np.linspace(0, 2 * np.pi, 100) 171 y = np.linspace(0, 2 * np.pi, 100).reshape( -1, 1) im = plt.imshow(f(x, y), cmap="hot", animated=True) def redraw(*args): global x, y x += np.pi / 26 y += np.pi / 20. im.set_array(f(x, y)) im.set_cmap('hot') return im, ani = animation.FuncAnimation(fig, redraw, interval=50, blit=True) plt.show() Один из кадров анимации показан на следующем рисунке. Напомним, что значит аргумент функции со звездочкой: redraw(*args). Если перед аргументом в определении функции указан символ * (звездочка), то функции можно передавать произвольное количество аргументов и они будут храниться в кортеже args. В приведенном коде функция redraw, вызываемая для рисования кадров, не использует свои аргументы. Но она меняет глобальные массивы x и y, которые в методе set_array используются для корректировки матрицы данных объекта im, созданного ранее функцией imshow. Обратите также внимание на то, как первоначально создается массив y. y = np.linspace(0, 2 * np.pi, 100).reshape(-1, 1) В методе reshape(-1,1) минус единица означает количество элементов по последней размерности. В результате создается двумерный массив с одним эдементом в каждой строке и количеством строк, равным количеству столбцов (элементов) одномерного массива. ■ Другой функцией, используемой для создания анимации, является функция ArtistAnimation . Она имеет сходные с функцией FuncAnimation аргументы, и отличие состоит только во втором. Он теперь является списком графических объектов, представляющих кадры анимации. Ниже приведен пример, использующий эту функцию. Пример. Последовательная анимация рисунков. В примере рисуются трехмерные объекты, показанные ниже. Из них создается список «кадров», которые затем по очереди отображаются в графической области. 172 Эти рисунки являются графиками поверхностей и пространственных кривых. Они строятся различными функциями, возвращающими значения в соответствующие переменные. Каждая из функций plot_surface возвращает объект класса mpl_toolkits.mplot3d.art3d.Poly3DCollection (коллекция 3D многоугольников), которые показанным в примере способом должны быть включены в список. Функция plot_wireframe возвращает объект класса Line3DCollection , а объекты, возвращаемые функциями plot и plot3D, являются пространственными кривыми. Все они включаются в анимационный список графических объектов похожим способом. Однако обратите внимание на то, как возвращаются объекты p6 и p7 (с запятой после имени объекта). Код программы приведен ниже. Большая его часть состоит из инструкций создания объектов рисунков, которые затем помещаются в список кадров, передаваемый для анимации в функцию ArtistAnimation. import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation from mpl_toolkits.mplot3d import Axes3D import copy fig = plt.figure(facecolor='white') ax=Axes3D(fig, xlim=( -1, 1), ylim=( -1, 1), zlim=( -1,1.0)) # поверхность 1 x = np.linspace( -1,1,31) y = np.linspace( -1,1,31) X,Y=np.meshgrid(x,y) Z=1-X**2-Y**2 ZZ=copy.deepcopy(Z) # запоминаем Z для 7-й поверхности p1=ax.plot_surface(X,Y,Z,rstride=1,cstride=1,color='r', alpha=0.8) # поверхность 2 Z=1-X**2-Y**2; Z[Z>0]=0 p2=ax.plot_surface(X,Y,Z,rstride=1,cstride=1,color='b', alpha=0.8) # поверхность 3 173 a=2.4; Z=np.sin(4*np.pi*(X**2+Y**2)/a**2) p3=ax.plot_surface(X,Y,Z,rstride=1,cstride=1,color='g', alpha=0.8) # поверхность 4 Z=1-X**2-Y**2; Z[Z<0]=0 Z=0.5-1.5*np.sqrt(Z) p4=ax.plot_surface(X,Y,Z,rstride=1,cstride=1,color='c', alpha=0.8) # поверхность 5 Z=1.3-np.abs(X-Y)-np.abs(X+Y) Z[Z<0]=0 Z=0.5-1.5*np.sqrt(Z) p5=ax.plot_wireframe(X,Y,Z,rstride=1,cstride=1,color='m', alpha=0.8, lw=2) # объект 6 (кривая) theta = np.linspace(-4 * np.pi, 4 * np.pi, 1000) z = np.linspace( -1, 1, 1000) r=np.abs(z) x = r * np.sin(6*theta) y = r * np.cos(6*theta) p6,=ax.plot(x, y, z,c='c',lw=3) # объект 7 (кривая в форме поверхности) # X,Y,Z первой поверхности ZZ[ZZ<0]=0 p7,=ax.plot3D(np.ravel(X),np.ravel(Y),np.ravel(ZZ),lw=2,c='brown') # список объектов для анимирования piclist=[(p1,),(p2,),(p3,),(p4,),(p5,),(p6,),(p7,)] anim = animation.ArtistAnimation(fig, piclist, interval=1000, blit=True, repeat_delay=100) plt.show() Обратите внимание на опцию repeat_delay=100. Без нее и с вкюченной опцией blit=True первый кадр (первая поверхность) на повторных циклах анимации не появляется (или незаметен). Если вы хотите, чтобы анимация не повторялась многократно, то используйте опцию repeat=False. ■ В следующем примере кадры анимации состоят из изображений матриц, которые «рисуются» функцией pcolor. Пример. Анимация изображений матриц. import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation import copy fig = plt.figure(facecolor='white') ax = plt.axes(xlim=( -25, 25), ylim=(-25, 25),aspect='equal') x = np.arange( -25, 26) 174 y = np.arange( -25, 26) X,Y=np.meshgrid(x,y) ZZ=np.sqrt(X**2+Y**2) imlist = [ ] for i in np.arange(36): Z=copy.deepcopy(ZZ) -i Z[np.logical_and(Z<0,Z> -2)]=(i-17)*2 Z[np.logical_and(Z<6,Z>4)]=(17 -i)*2 imlist.append((plt.pcolor(x, y, Z , norm=plt.Normalize( -35, 35)),)) anim = animation.ArtistAnimation(fig, imlist, interval=50, repeat_delay=100, blit=True) plt.show() Несколько кадров, созданных этой программой, приведены на следующем рисунке. Здесь список imlist составляется из графических объектов, возвращаемых функцией pcolor. Она вызывается в одном из двух форматов: pcolor(C[,...]) pcolor(X, Y, C[, ...]) Здесь С – массив значений, которые представляют цвета. Если используются массивы X и Y, то они определяют (x, y) координаты закрашиваемых четырехугольников, вершины которых расположены в точках (X[i,j],Y[i,j]), (X[i,j+1],Y[i,j+1]), (X[i+1,j],Y[i+1,j]), (X[i+1,j+1], Y[i+1,j+1]). При этом индекс столбца ссответствует x координате, а индекс строки – y координате. Размеры массивов X и Y должны быть на единицу больше соответствующих размеров матрицы C. Однако, если размеры одинаковы, то последняя строка и столбец матрицы C игнорируются. В качестве x и y массивов допустимо использовать одномерные массивы. Аргумент norm=plt.Normalize(-35,35) функции pcolor используется для нормировки данных (приведению их к диапазону [0, 1]). 4.4 Графические функции модуля mpmath В набор научных библиотек входит пакет mpmath, который позволяет выполнять вычисления стандартных математических функций (элементарных и специальных) с произвольной точностью. В нашем пособии мы этой темы не будем касаться. Но в этот модуль включено несколько функций, которые умеют строить графики. Их особенность состоит в том, что для построения графиков функций требуются только их аналитические представления. 175 Функция mpmath.plot(...) строит график одной или нескольких функций одной вещественной переменной. При этом используется идентификатор имени функции (без аргументов). from mpmath import * plot([cos,sin],[0,4*pi]) # следующий рисунок слева f=lambda x: sin(x)/x plot(f, [-24, 24]) # предыдущий рисунок справа Имеются опции, управляющие пределами изменения значений по осям, а также опция singularities, которая указывает особые точки графика. f=lambda x: tan(x)*x plot(f,xlim=[ -6,6],ylim=[-10,10],singularities=[ -3*pi/2,-pi/2, pi/2,3*pi/2]) Функция mpmath.cplot(f,re=[-5,5],im=[-5,5],points=2000, color=None,...) строит график комплекснозначной функций f над прямоугольной областью комплексной плоскости, задаваемой парой интервалов re и im. По умолчанию аргумент функции представляется цветом в hue палитре, а модуль значения представляется яркостью. cplot(lambda z: z, [-2,2], [-2,2]) # следующий рисунок слева cplot(lambda z: z**3, [-2,2], [-2,2]) # предыдущий рисунок справа Вот как выглядит график функции Жуковского z z z f 1 2 1 176 cplot(lambda z: (z+1/z)/2, [-2,2], [-2,2]) # следующий рисунок слева На предыдущем рисунке справа показан график функции 1 2 w w z обратной к функции Жуковского. cplot(lambda w: w+sqrt(w**2-1),[-2,2],[-2,2],points=20000) Функция mpmath.splot(f,u=[-5,5],v=[-5,5],points=100, keep_aspect=True, wireframe=False,...) строит график (поверхность) функции двух переменных. f=lambda x, y: cos(x**2 -y**2)*exp(-(x**2+y**2)/16) splot(f, [-3,3], [-3,3]) # следующий рисунок слева f=lambda x, y: abs(x -y)+abs(x+y) splot(f, [-2,2],[-2,2],points=200,wireframe=True) # пред. рисунок справа Если f возвращает список из трех функций, то функция splot() рассматривает их как компоненты параметрического уравнения поверхности. r=0.5 f=lambda ph,th: [(1+r*cos(ph))*cos(th),(1+r*cos(ph))*sin(th),r*sin(ph)] splot(f,[0,2*3.1415], [0,2*3.1415]) # следующий рисунок слева 177 R,H=1,3 # радиус и высота конуса fx=lambda t,tau: R*cos(tau)*(abs(t) -2*abs(t-1)+abs(t-2))/2 fy=lambda t,tau: R*s in(tau)*(abs(t) -2*abs(t-1)+abs(t-2))/2 fz=lambda t,tau: H*(1+abs(t -1)-abs(t-2))/2 f=lambda t,tau: [fx(t,tau),fy(t,tau),fz(t,tau)] splot(f,u=[0,2],v=[0,2*3.141592],wireframe=True,points=200) В этом примере мы по параметрическим уравнениям нарисовали конус, который показан на предыдущем рисунке справа. 5. Символьные вычисления 5.1 Основы символьных вычислений Python может выполнять арифметические операции не только как калькулятор. После импортирования функций модуля SymPy, он может выполнять символьные преобразования. >>> from sympy import * Поскольку переменные в Python не имеют никакого значения до тех пор, пока им не будет присвоено значение, то просто использовать их в символьных вычислениях нельзя. Поэтому перед использованием переменные и функции должны быть объявлены как символьные. Это отличает Python SymPy от таких программ как Mathematica или Maple, где переменные, если им не присвоено никаких значений, автоматически рассматриваются как символьные. Объявление символьных переменных выполняется функцией symbols. Например >>> x,y,a,b = symbols('x y a b') Здесь мы создали четыре символьных переменных. Предыдущие значения переменных затираются. Символьные выражения конструируются из символьных переменных >>>f=a**3*x + 3*a**2*x**2/2 + a*x**3 + x**4/4 Переменная f автоматически становится символьной. Имеется другая функция var(), которая также создает символьные переменные, помещая их в глобальное пространство имен. >>> var('u,v') (u, v) Разница между symbols() и var() в том, первая функция возвращает ссылку на символьный объект, которую, чтобы использовать, нужно присвоить какой – либо переменной. Вторая, без присваивания создает символьную переменную. Например, после выполнения команды >>>var('Hello') символьная переменная с именем Hello будет уже существовать. >>> var('Hello') Hello >>> type(Hello) При создании символьных выражений следует аккуратно обращаться с числовыми константами. Иногда они могут быть проинтерпретированы как 178 числовые константы Python, а не SymPy. Если надо быть уверенным, что константа является символьной, поместите ее внутрь функции S(). >>> expr = x**2 + sin(y) + S(1)/2; expr x**2 + sin(y) + 1/2 Сравните >>> type(1) >>> type(S(1)) # символьная константа единица Функция sympy.S(expr[,...]) конвертирует произвольное выражение expr в тип, который можно использовать в выражениях SymPy. Например, она преобразует целые, рациональные, вещественные числа Python в аналогичные типы SymPy. Ее аргумент expr может быть любым объектом SymPy, объектом стандартного Python типа (int, long, float, Decimal), строкой, булевской переменной, списком, множеством или кортежем из элементов приведенных типов. >>> type(S(x**2)) sympy.core.power.Pow >>> S(x**2).subs(x,3) 9 Также функцию S() можно использовать и для объявления символьных переменных. >>> p = S('p') >>> type(p) sympy.core.symbol.Symbol Фактически инструкция S() является синонимом вызова функции sympify(), которая часто используется для преобразования строки в символьное выражение. Но, как сказано выше, ее (функцию S()) можно использовать также для преобразования констант Python в символьные константы. >>> sympify(0.5) 0.500000000000000 >>> type(_) # тип последнего результата Разница между константой Python и символьной состоит в том, что символьная постоянная может быть вычислена с любой степенью точности. Для этого можно использовать метод симв_перем.n(кол_цифр). Например, следующая инструкция вычисляет переменную z1 с процессорной точностью (примерно 15 значащих цифр). >>> z1=1/3; z1 0.3333333333333333 Символьная константа может быть вычислена с произвольной точностью. >>> z2=S(1)/3; z2 1/3 >>> z2.n(30) 0.333333333333333333333333333333 179 Сравните также следующие два выражения. >>> x**(1/2) x**0.5 >>> x**(S(1)/2) sqrt(x) Можно объявлять символьные переменные с индексом. >>> x=symbols('x:5'); x # диапазон индексов от 0 до 4 (x0, x1, x2, x3, x4) #x[0],x[1],... явл. симв. переменными x0,x1,... >>> x=symbols('x5:10'); x # диапазон индексов от 5 до 9 (x5, x6, x7, x8, x9) Можно создать набор идущих подряд «однобуквенных» символьных имен. >>> symbols('a:d') (a, b, c, d) При создании символьным переменным можно назначить тип, например, целый. >>> k, m, n = symbols('k m n', integer=True) Иногда без подобных дополнительных предположений почти очевидные символьные преобразования не работают. >>> sqrt(x**2) # это x, если x≥0 sqrt(x**2) Функция Symbol(„имя‟,...) (на самом деле это конструктор класса) позволяет наложить на символьную переменную какое–нибудь условие. >>> x = Symbol('x', positive=True) >>> sqrt(x**2) x Функция var() также имеет опции наложения ограничений на символьные переменные. >>> var('x y', positive = True, integer = True) Если вы предполагаете в текущей сессии использовать символьную математику постоянно, то можно импортировать общепринятые символьные имена из модуля sympy.abc. >>> import sympy.abc >>> dir(sympy.abc) ['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', ..., 'a','alpha','b','beta','c', 'chi', 'd','delta', 'division', 'e', 'epsilon', 'eta', 'exec_', 'f', 'g', 'gamma', 'greeks', 'h', 'i', 'iota', 'j', 'k', 'kappa','l','lamda', 'm', 'mu', 'n', 'nu', 'o', 'omega','omicron','p', 'phi', 'pi', 'print_function', 'psi', 'q', 'r', 'rho', 's', 'sigma', 'string', 'symbols', 't', 'tau', 'theta', 'u', 'upsilon', 'v', 'w', 'x', 'xi', 'y', 'z', 'zeta'] Если импортировать содержимое модуля sympy.abc командой >>> from sympy.abc import * то все имена, перечисленные выше (мы несколько имен пропустили, заменив их символом многоточия) можно использовать в символьных выражениях. 180 Однако мы не рекомендуем это делать, поскольку будут перегружены имена некоторых символьных констант, например, таких как pi или E (основание натурального логарифма). Напомним, что удалить имя переменной из пространства имен можно командой del имя1, имя2,... >>> type(x) >>> del x,y >>> x Ошибка! Если вы выполняли инструкцию « from sympy.abc import * », то повторно загрузите модуль sympy. >>> from sympy import * Этим вы восстановите значения стандартных констант pi и E, а также имена некоторых функций. При записи символьного выражения может автоматически выполняться его упрощение. >>> a,b,c,d,x,y,z,u,v,w = symbols('a b c d x y z u v w') >>> x - z + 23 -z- 14 + sin(pi)+2*z x + 9 Для вычисления символьного выражения при некоторых конкретных значениях переменных используется метод subs(...). >>> f=a**3*x + 3*a**2*x**2/2 + a*x**3 + x**4/4 >>> f.subs(a,1) x**4/4 + x**3 + 3*x**2/2 + x Здесь в выражение f вместо переменной a была подставлена единица. Если метод subs принимает два аргумента, то они интерпретируются как subs(old,new) , т.е. старый идентификатор old заменяется новым new. Аргумент метода subs() может быть последовательностью, которая должна содержать пары (old,new). В следующей команде в том же символьном выражении f выполнена двойная подстановка a=1, x=2. >>> f.subs([(a,1),(x,2)]) 20 >>> f=x*(1+log(x)) >>> f.subs(x, pi).evalf() 6.73786765331895 Подставляемый аргумент может быть словарем, чьи пары ключ:значение должны соответствовать старым и новым значениям old:new. >>> f.subs({x:1}) a**3 + 3*a**2/2 + a + 1/4 Вместо символьной переменной подставлять можно символьное выражение. >>> expr=x**2+2*x+1 >>> expr.subs(x,1/x) 1 + 2/x + x**(-2) 181 Обратим ваше внимание на следующую особенность работы с переменными (символьными и обычными переменными Python). Выполним следующий код >>> x = symbols('x') >>> expr=x+1 >>> x=2 >>> print(expr) x + 1 Присваивание символьной переменной x значения 2 (этим меняем также тип x) не повлияло на значение символьного выражения expr. Здесь действует правило: если переменная изменилась, то созданное ранее выражение, содержащее эту переменную, не пересчитывается автоматически. Это правило срабатывает и для обычных переменных Python. Например, >>> x='Hello' >>> expr=x+' world' >>> expr 'Helloworld' >>> x='AAA' >>> expr 'Helloworld' Символьные вычисления также имеют отношение к выражениям, в которые входят одни числа. Например, SymPy может проводить вычисления с дробями и, приводя их к общему знаменателю, и получать точный ответ. >>> S(1)/3+S(2)/5 11/15 Сравните. >>> 1/3+2/5 0.7333333333333334 Для создания рациональных дробей можно использовать функцию Rational(числитель,знаменатель). С рациональными числами математические операции выполняются без десятичного округления. >>> z=Rational(1, 3)+Rational(2, 5); z 11/15 Для создания рациональной дроби можно также использовать функцию Integer(...). Например, >>> Integer(1)/Integer(3) 1/3 >>> 1/3 0.3333333333333333 >>> z=Integer(1)/Integer(3)+Rational(2, 5); z 11/15 Основной принцип символьных вычислений – не делать никаких приближений, кроме затребованных. Так выражение >>> sqrt(12) # функция sqrt() модуля sympy будет преобразовано к виду 2*sqrt(3) 182 Python преобразует выражение, но оставляет в записи ответа квадратный корень и не выполняет никаких округлений. Однако, если хотя бы одно из чисел символьного выражения задано с десятичной точкой, то и результат будет приближенным >>> sqrt(12.0) 3.46410161513775 У любого символьного объекта имеется метод evalf(...)( |