NumProblems

Материал из theor
Перейти к: навигация, поиск

[Пример Адамара (некорректно поставленная задача)|http://ru.wikipedia.org/wiki/%D0%9F%D1%80%D0%B8%D0%BC%D0%B5%D1%80_%D0%90%D0%B4%D0%B0%D0%BC%D0%B0%D1%80%D0%B0]

Python

Маятник Капицы

Формулировка задачи

Рассчитать численно движение маятника с подвижным основанием, проверить случай устойчивого "перевёрнутого" положения.

  • Общее описание

→ Система представляет собой точечный груз, соединённый с основанием невесомым абсолютно твёрдым стержнем. При этом, основание маятника совершает периодическое колебательное движение.
→ Рассмотрим случай, когда основание маятника колеблется вдоль вертикальной оси по закону z=a*sin(wt).

  • Параметры и переменные

g - ускорение свободного падения
m a - амплитуда колебаний основания
l - длина стержня
w - частота колебаний основания
t - время
φ - угол поворота стержня, отсчитываемый от направления g

  • Используемые библиотеки

пакет NumPy - для работы с массивами.
модуль Math - для вычисления тригонометрических функций.
интерфейс PyLab (библиотека Matplotlib) - для рисования графика.

Аналитическая часть решения

  • Запишем лагранжиан системы <math>\begin{matrix} L = T + U \end{matrix}</math>, где
<math>\begin{matrix}

T = \frac{m l^2 }{2} \dot \varphi^2 + m a l w ~\sin(w t) \sin(\varphi)~\dot\varphi + \frac{m a^2 w^2}{2} \sin^2(w t)

\end{matrix}</math>
- кинетическая энергия маятника,
<math>\begin{matrix}

U = -m g (l cos(\varphi) + a cos(w t))

\end{matrix}</math>
- его потенциальная энергия.
  • Отсюда, записав уравнение Лагранжа, получаем:
<math>\begin{matrix}

\ddot \varphi = - (a~w^2\cos w t~ +g) \frac{\sin \varphi}{l}\;,

\end{matrix}</math>

Это уравнение и будем решать численно. В результате получим зависимость угла φ от времени.

Текст программы

<source lang = 'python'>

 import  math 
 import numpy
 import pylab
 g=10
 l=1
 w=250
 a=0.1
 #Начальное положение маятника:
 fi=3.1415/2+1
 #Начальная скорость маятника:
 fit=0
 t=0
 dt=0.001
 n=0
 nmax=1000
 # Готовим векторы для графиков:
 tt=import numpy.zeros(nmax+1)
 fifi=import numpy.zeros(nmax+1)
 zz=import numpy.zeros(nmax+1)
 xx=import numpy.zeros(nmax+1)
 while n<nmax+1:
 ::	F=-1*math.sin(fi)*(1/l)*(g+a*w*w*math.sin(w*t))
 ::	fit=fit+F*dt
 ::	fi=fi+fit*dt
 ::	zz[n]=a*math.cos(w*t)-l*math.cos(fi)
 ::	xx[n]=l*math.sin(fi)
 ::	fifi[n]=fi
 ::	tt[n]=t
 ::	t=t+dt
 ::	n=n+1
 
 pylab.plot(xx,zz)
 pylab.show()

</source>
Шаг по времени dt должен удовлетворять условию: <math>\begin{matrix} wdt<<2\pi \end{matrix}</math>, где w - максимальная частота процессов в системе.

Результаты деятельности

Траектория маятника в плоскости x,y

Рассчитав через полученную φ(t) зависимость координаты z от времени (как более наглядную величину), забиваем в программу разные условия. При достаточно больших значениях частоты колебаний основания и амплитуде a << l, приходим к задаче, сформулированной П.Л. Капицей для иллюстрации придуманного им метода эффективного потенциала.
Если теперь начальное значение φ положить близким к π, то получим колебания относительно "перевёрнутого" положения равновесия.
Следует отметить, что мы не накладывали никаких ограничений на амплитуду и частоту колебаний основания.

Дополнительные задания

  • Рассчитать аналогичную задачу для случайных колебаний основания (случайная сила). Будет ли в этом случае существовать устойчивое положение равновесия в "перевёрнутом" состоянии?

Квантовые часы

Описание - Квантовые часы

Текст программы (метод Дюфорта-Франкеля)

<source lang="python"> from pylab import * from numpy import * from matplotlib import animation dx = 0.2 dt = 0.02 mult = round(0.3 / dt) g = 1j*dt/dx**2 xmax = 40.0 sp_x = arange(-xmax,xmax,dx) sp_t = arange(0.0,100.0,dt) nx = size(sp_x) b = 1.0 a = 1.0 v0 = -1.0 p0 = 3.0

  1. Начальные условия
  1. Гауссов пакет с импульсом
  2. init = array(map(lambda x:exp(-(x+10)**2/(2*(a)**2)+1j*p0*x),sp_x))
  1. Покоящаяся частица в правой яме

init = array(map(lambda x: exp(-(x-(b+a/2))**2.0/(2*a)**2),sp_x), dtype=dtype('complex'))

  1. Нормировка

init = init/sqrt(dot(init,conj(init)))


def v(x): return v0/cosh((x-(b+a/2.0))/a)**2 + v0/cosh((x+(b+a/2.0))/a)**2

V = array(map(v,sp_x),dtype=dtype('complex')) u = zeros((size(sp_t),size(sp_x)), dtype=dtype('complex')) u[0,:] = init

  1. first step - explicit scheme

for x in range(size(sp_x)): if (x!=0) and (x != size(sp_x)-1): u[1,x]=u[1-1,x]+(1j)*dt*(u[1-1,x+1] - 2*u[1-1,x] + u[1-1,x-1])/dx**2 - 1j*dt*V[x]*u[0,x] if x==0: u[1,x]=u[1-1,x]+(1j)*dt*(u[1-1,x+1] - 2*u[1-1,x])/dx**2 - 1j*dt*V[x]*u[0,x] if x==size(nx)-1: u[1,x]=u[1-1,x]+dt*(1j)*( -2*u[1-1,x] +u[1-1,x-1])/dx**2 - 1j*dt*V[x]*u[0,x]

  1. Duforte-Frankel

import progressbar pbar = progressbar.ProgressBar(maxval = size(sp_t)).start() import sys for t in range(2,size(sp_t)): pbar.update(t)

for x in range(size(sp_x)): if (x!=0) and (x!=size(sp_x)-1): temp = (2j*dt/(dx**2.0))*(u[t-1,x+1] + u[t-1, x-1] - u[t-2, x]) temp += -2j*dt*V[x]*u[t-1,x] + u[t-2,x] u[t,x] = temp/(1.0+2j*dt/dx**2.0) if (x==0): temp = (2j*dt/dx**2.0)*(u[t-1,x+1] - u[t-2, x]) temp += -2j*dt*V[x]*u[t-1,x] + u[t-2,x] u[t,x] = temp/(1.0+2j*dt/dx**2.0) if (x==size(sp_x) - 1): temp = (2j*dt/dx**2.0)*(u[t-1, x-1] - u[t-2, x]) temp += -2j*dt*V[x]*u[t-1,x] + u[t-2,x] u[t,x] = temp/(1.0+2j*dt/dx**2.0)


fig = plt.figure() ax = fig.add_subplot(111) line, = ax.plot(sp_x, map(lambda x: abs(x)**2,init)) ax.plot(sp_x,V)

  1. ax.plot(sp_x, u[1,:])
  2. show()

def animate(i):

   line.set_ydata(map(lambda x: abs(x)**2,u[mult*i,:]))
   return line

interval = 1

ani = animation.FuncAnimation(fig,animate, np.arange(1,size(sp_t)/mult), interval=interval) plt.show() </source>

Результаты работы

Сами квантовые часы
: thubm

Дополнительные задания

  • Посмотреть на туннелирование пакета через барьер
  • ...TODO...

    Уравнение теплопроводности

    Описание - Задача теплопроводности

    Программа на Python

    <source lang = 'python'>

    1. -*- coding: utf-8 -*-

    from time import * from pylab import * from matplotlib import pyplot as plt from matplotlib import animation

    class Grid:

       def __init__(self, init, bc_left, bc_right, nt = 50, nx = 50, xmin = 0.0, 
                    xmax = 1.0, tmin = 0.0, tmax = 5.0):
           self.nx = nx
           self.nt = nt
           self.xmin = xmin
           self.xmax = xmax
           self.tmin = tmin
           self.tmax = tmax
           self.x, self.dx = linspace(xmin, xmax, nx, retstep=1)
           self.t, self.dt = linspace(tmin, tmax, nt, retstep=1)
           if (self.dt > self.dx**2):
               print u'Явная схема неустойчива!'
           self.u = zeros((nt, nx), dtype='f')
           self.u[0, :] = map(init, self.x)
           self.u[:, 0] = map(bc_left, self.t)
           print self.u[:,0]
           self.u[:, self.nx - 1] = map(bc_right, self.t)
           
    

    class Solver:

       def __init__(self, grid, stepper = 'explicit'):
           self.grid = grid
           self.setTimeStepper(stepper)
           
       def slowTimeStep(self):
           dt = self.grid.dt
           dx = self.grid.dx
           u = self.grid.u
           nt, nx = self.grid.u.shape
           t = time()
           for t in range(1,nx-1):
               for x in range(1, nx-1):
                   u[t,x] = u[t-1,x] + dt*(u[t-1,x+1] - 2*u[t-1,x] +u[t-1,x-1])/dx**2
           return time() - t
           
       def explicitTimeStep(self):                    
           dt = self.grid.dt
           dx = self.grid.dx
           u = self.grid.u
           t = time()
           u[1:-1, 1:-1] = u[0:-2, 1:-1] + dt*(u[0:-2,2:] - 2*u[0:-2,1:-1] + u[0:-2,0:-2])/dx**2
           return time() - t
           
       def implicitTimeStep(self):
           return 0;
       
       def setTimeStepper(self, stepper = 'explicit'):
           if stepper == 'explicit_slow':
               self.timeStep = self.slowTimeStep
           if stepper == 'explicit':
               self.timeStep = self.explicitTimeStep
           if stepper == 'implicit':
               self.timeStep = self.implicitTimeStep
       
       def Solve(self):
           t = time()
           for i in xrange(1, self.grid.nt):
               print i
               self.timeStep()
           print time()-t
           return self.grid.u
       
       def ani(self, i):
           self.line.set_ydata(map(lambda x: x,self.grid.u[i,:]))
           return self.line
       
       def animate(self, interval = 30):
           u = self.grid.u
           fig = plt.figure()
           ax = fig.add_subplot(111)
           self.line, = ax.plot(self.grid.x, map(lambda x: abs(x)**2,u[0,:]))
           ani = animation.FuncAnimation(fig, self.ani,self.grid.nt, interval=interval)
           plt.show()
           
           
           
           
    

    g = Grid(lambda x: exp(-x**2), lambda x: 0.0, lambda x: 0.0, xmin = -1.0, nx = 30, nt = 4000) S = Solver(g, stepper='explicit_slow') U = S.Solve() S.animate

    </source> В данной программе использованы два способа присваивания - в цикле и при помощи срезов NumPy (сравните скорость!).Неявную схему предлагается скопипастить самостоятельно по типу: <source lang='python'>

           M = zeros((nx,nx))#матрица коэффициентов системы
    

    for j in range(1, nx-1): M[j,j+1] = g M[j,j] = -(1.0+2.0*g) M[j,j-1] = g res = linalg.solve(M,y)#решение её return res </source>


    Потенциал

    <source lang = 'python'>

    1. -*- coding: utf-8 -*-

    """ Created on Sat Jun 01 10:49:09 2013

    """

    import numpy from mpl_toolkits.mplot3d.axes3d import Axes3D import matplotlib.pyplot as plt from matplotlib import cm from mpl_toolkits.mplot3d.axes3d import get_test_data

    class Grid:

      def __init__(self, nx=10, ny=10, xmin=0.0, xmax=1.0,
                    ymin=0.0, ymax=1.0):
           self.xmin, self.xmax, self.ymin, self.ymax = xmin, xmax, ymin, ymax
           self.dx = float(xmax-xmin)/(nx-1)
           self.nx = nx
           self.ny = ny
           self.dy = float(ymax-ymin)/(ny-1)
           self.u = numpy.zeros((nx, ny), 'd')
           self.u[nx/2,ny/2] = -1.0
           # used to compute the change in solution in some of the methods.
           self.old_u = self.u.copy()
    
       def setBCFunc(self, func):
           xmin, ymin = self.xmin, self.ymin
           xmax, ymax = self.xmax, self.ymax
           x = numpy.arange(xmin, xmax + self.dx*0.5, self.dx)
           y = numpy.arange(ymin, ymax + self.dy*0.5, self.dy)
           self.u[0 ,:] = func(xmin,y)
           self.u[-1,:] = func(xmax,y)
           self.u[:, 0] = func(x,ymin)
           self.u[:,-1] = func(x,ymax)
           
    
       def computeError(self):
         
           v = (self.u - self.old_u).flat
           return numpy.sqrt(numpy.dot(v,v))
    


    class LaplaceSolver:

       def __init__(self, grid, stepper='numeric'):
           self.grid = grid
           self.setTimeStepper(stepper)
    
       def slowTimeStep(self, dt=0.0):
           g = self.grid
           nx, ny = g.u.shape
           dx2, dy2 = g.dx**2, g.dy**2
           dnr_inv = 0.5/(dx2 + dy2)
           u = g.u
    
           err = 0.0
           for i in range(1, nx-1):
               for j in range(1, ny-1):
                   tmp = u[i,j]
                   u[i,j] = ((u[i-1, j] + u[i+1, j])*dy2 +
                             (u[i, j-1] + u[i, j+1])*dx2)*dnr_inv
                   diff = u[i,j] - tmp
                   err += diff*diff
    
           return numpy.sqrt(err)
       
       def numericTimeStep(self, dt=0.0):
          
           g = self.grid
           dx2, dy2 = g.dx**2, g.dy**2
           dnr_inv = 0.5/(dx2 + dy2)
           u = g.u
           g.old_u = u.copy() 
    
           # The actual iteration
           u[1:-1, 1:-1] = ((u[0:-2, 1:-1] + u[2:, 1:-1])*dy2 +
                            (u[1:-1,0:-2] + u[1:-1, 2:])*dx2)*dnr_inv
    
           return g.computeError()
    
       def setTimeStepper(self, stepper='numeric'):
           
           if stepper == 'slow':
               self.timeStep = self.slowTimeStep
           # ...
           else:
               self.timeStep = self.numericTimeStep
    
       def solve(self, n_iter=0, eps=1.0e-16):
           err = self.timeStep()
           count = 1
    
           while err > eps:
               if n_iter and count >= n_iter:
                   return err
               err = self.timeStep()
               count = count + 1
    
           return count
    

    g = Grid(nx = 100, ny = 100) S = LaplaceSolver(g) print S.solve(100)

    fig = plt.figure(figsize=plt.figaspect(0.5)) ax = fig.add_subplot(1, 2, 1, projection='3d') X = numpy.linspace(S.grid.xmin, S.grid.xmax, S.grid.nx ) Y = numpy.linspace(S.grid.ymin, S.grid.ymax, S.grid.ny ) X,Y = numpy.meshgrid(X,Y) Z = S.grid.u surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,

           linewidth=0, antialiased = False)
           
    

    plt.show() </source>

    Решение уравнения Кортевега-де Фриза

    Текст программы

    Файл:Решение уравнения Кортевега-де Фриза на Python.pdf