Как решать дифференциальные уравнения в матлабе
Перейти к содержимому

Как решать дифференциальные уравнения в матлабе

  • автор:

Решение ОДУ в Matlab

Решение ОДУ в Matlab

Доброго времени суток! Сегодня мы поговорим о решении ОДУ (обыкновенных дифференциальных уравнений) в Matlab. Перед тем как мы начнём обсуждать данную тему, советую вам ознакомиться с темой: Численное дифференцирование в Matlab, чтобы лучше понимать теоретическую составляющую решения ОДУ.

Обыкновенные дифференциальные уравнения

  • Задача Коши
  • Краевая задача
  • Задача на собственные значения

Кратко расскажу о их сути:

Задача Коши предполагает дополнительные условия в виде значения функции в определённой точке.
Краевая задача подразумевает поиск решения на заданном отрезке с краевыми (граничными) условиями в концах интервала или на границе области.
Задача на собственные значения — помимо искомых функций и их производных, в уравнение входят дополнительное несколько неизвестных параметров, которые являются собственными значениями.

Методы решения дифференциальных уравнений

Решение ОДУ в Matlab и не только, в первую очередь, сводится к выбору порядка численного метода решения. Порядок численного метода не связан с порядком дифференциального уравнения. Высокий порядок у численного метода означает его скорость сходимости.

В случае большого интервала, с помощью алгоритмов с низким порядком сжимают интервал с решениями и находят приблизительные корни, а затем уже уточняют корни с помощью методов с высоким порядком.

Решение обыкновенных дифференциальных уравнений в Matlab можно реализовать «своими ручками», прописав алгоритм по разным схемам. Но также в Matlab есть встроенные функции, выполняющие все стандартные задачи.

Метод Рунге-Кутта первого порядка

Методы Рунге-Кутта представляют собой разложения в ряд Тейлора и от количества использованных элементов ряда зависит порядок этого метода. Следовательно, помимо Рунге-Кутта первого порядка, вы сможете увидеть методы других порядков. Иногда их называют другими именами.

Например, Метод Рунге-Кутта первого порядка, также известен как Метод Эйлера или Метод ломаных. Информацию о его математическом и графическом представлении советую поискать в гугл. Мы же поговорим о том, как Метод Рунге-Кутта первого порядка реализуется в Matlab для решения ОДУ. Например:

Решить и привести график ошибки уравнения y’ = y*x методом Рунге-Кутта первого порядка (Методом Эйлера, Методом ломаных).

Setka=10:10:10000; for k=1:length(Setka) %определяем параметры сетки N=Setka(k); h=1.0/(N-1); %задаем начальное условие y(1)=1; %применяем алгоритм метода ломаных/Эйлера for n=1:(N-1) y(n+1)=y(n)+ h*((n-1)*h*y(n)); %где (n-1)*h -> x end %вычисляем величину M(1) в оценке %погрешности численного решения M(k)=abs(y(N)-exp(0.5))/h; step(k)=h; end %рисуем зависимость величины M(1) от шага plot(step,M); 

Погрешность метода 1 порядка

На данном графике показана зависимость величины ошибки от шага.

Метод Рунге-Кутта второго порядка

Также известен как Метод Эйлера-Коши. Как видите, во второй части уравнения происходит обращения к следующему шагу. Но как тогда быть, если нам ещё не известен следующий шаг? Всё просто. Метод Рунге-Кутта второго порядка — это всё тот же метод первого порядка, однако, на половине шага происходит нахождение «первичного» решения, а затем происходит его уточнение. Это позволяет поднять порядок скорости сходимости до двух.

Решить и привести график ошибки уравнения u’ = u*x методом Рунге-Кутта второго порядка.

f=@(x,u)x*u; %задаем набор сеток Setka=10:50:10000; %организуем цикл расчетов с разными сетками for k=1:length(Setka) N=Setka(k); h=1.0/(N-1); %задаем начальное условие y(1)=1; %вычисляем приближенные значения решения %дифференциального уравнения for n=1:(N-1) x(n)=(n-1)*h; y(n+1)=y(n)+h*(0.75*f(x(n),y(n))+... 0.25*f(x(n)+h/0.5,y(n)+... (h/0.5)*f(x(n),y(n)))); end %находим константу M в оценке погрешности %приближенного решения |y(N)-u(N)| %она не должна зависеть от h M(k)=abs(y(N)-exp(0.5))/h^2; step(k)=h; end %рисуем зависимость константы M от шага сетки h plot(step,M); 

Погешность Рунге-Кутта второго порядка

По сравнению с Рунге-Куттом первого порядка изначальная ошибка уже гораздо меньше.

Мы не будем говорить о третьем порядке, потому что задачи на третий порядок встречаются редко, но если будет необходимо, пишите в комментариях, выложу.

Метод Рунге-Кутта четвёртого порядка

Метод Рунге-Кутта четвёртого порядка считается самым распространённым. Тем не менее, работает он аналогично второму и третьему порядку.

Решить и привести график ошибки уравнения u’ = u*x методом Рунге-Кутта четвёртого порядка.

f=@(x,u)x*u; %задаем набор сеток Setka=10:50:10000; %организуем цикл расчетов с разными сетками for k=1:length(Setka) N=Setka(k); h=1.0/(N-1); %задаем начальное условие y(1)=1; %вычисляем приближенные решения %дифференциального уравнения for n=1:(N-1) x(n)=(n-1)*h; k1=f(x(n),y(n)); k2=f(x(n)+0.5*h,y(n)+0.5*h*k1); k3=f(x(n)+0.5*h,y(n)+0.5*h*k2); k4=f(x(n)+h,y(n)+h*k3); y(n+1)=y(n)+(h/6)*(k1+2*k2+2*k3+k4); end %находим константу M в оценке погрешности %приближенного решения |y(N)-u(x(N))| %она не должна зависеть от h M(k)=abs(y(N)-exp(0.5))/h^4; step(k)=h; end %рисуем зависимость константы M от шага сетки h loglog(step,M); 

Зависимость погрешности от шага

Как видите, на последней картинке размерность ошибки на столько мала, что пришлось воспользоваться loglog() для лучшей видимости.

Решение ОДУ в Matlab стандартными средствами

Стоит отметить, что мы с вами разобрали только один самый известный метод решения ОДУ с разными порядками. Однако, методов очень много.

Для решения дифференциальных уравнений и систем в MATLAB предусмотрены следующие функции:

ode45 (f, interval, X0, [options])
ode23 (f, interval, X0, [options])
ode113 (f, interval, X0, [options])
ode15s (f, interval, X0, [options])
ode23s (f, interval, X0, [options])
ode23t (f, interval, X0, [options])
ode23tb (f, interval, X0, [options])

Входными параметрами этих функций являются:

  • f — вектор-функция для вычисления правой части уравнения системы уравнений;
  • interval — массив из двух чисел, определяющий интервал интегрирования дифференциального уравнения или системы;
  • Х0 — вектор начальных условий системы дифференциальных уравнений;
  • options — параметры управления ходом решения дифференциального уравнения или системы.

  • массив Т — координаты узлов сетки, в которых ищется решение;
  • матрицу X, i-й столбец которой является значением вектор-функции решения в узле Тi.

В функции ode45 реализован метод Рунге-Кутта 4-5 порядка точности, в функции ode23 также реализован метод Рунге-Кутта, но 2-3 порядка, а функция ode113 реализует метод Адамса.

Для решения жёстких систем предназначены функция ode15s, в которой реализован метод Гира, и функция ode23s, реализующая метод Розенброка. Для получения более точного решения жёсткой системы лучше использовать функцию ode15s. Для решения системы с небольшим числом жёсткости можно использовать функцию ode23t, а для грубой оценки подобных систем служит функция ode23tb.

Символьное решение обыкновенных дифференциальных уравнений произвольного порядка осуществляет функция dsolve r = dsolve(‘eq1,eq2,…’, ‘cond1,cond2,…‘, ‘v’)
Пример использования:

% в отдельном М-файле с именем pr10 function f=pr10(x,y) f=sin(x+y)+(3/2)*(x-y); end) %в основном М-файле ode113(@pr10,[0 20],0) %Метод Адамса: @pr10 – ссылка на М-функцию, % [0 20]- интервал, 0 - условие: y(0)=0 
%пример (u'=u^2) expl=dsolve('Du=u^2','x') 

На этом мы закончим. Если остались вопросы, задавайте их в комментариях. Также вы можете скачать исходники чтобы лучше понять тему: «Решение ОДУ в Matlab».

Поделиться ссылкой:

Документация

Решение дифференциальных уравнений с частными производными

В partial differential equation (УЧП) функция, решаемая для, зависит от нескольких переменных, и дифференциальное уравнение может включать частные производные, взятые относительно каждой из переменных. Дифференциальные уравнения с частными производными полезны для моделирования волн, теплового потока, жидкой дисперсии и других явлений с пространственным поведением, которое изменяется в зависимости от времени.

Какие типы УЧП можно решить с MATLAB ?

MATLAB ® Решатель УЧП pdepe решает начальные краевые задачи для систем УЧП в одной пространственной переменной x и время t . Можно думать о них как об ОДУ одной переменной, которые также изменяются относительно времени.

pdepe использует неофициальную классификацию для 1D уравнений, которые она решает:

  • Уравнениями с производной времени является parabolic. Примером является уравнение тепла ∂ u ∂ t = ∂ 2 u ∂ x 2 .
  • Уравнениями без производной времени является elliptic. Примером является уравнение Лапласа ∂ 2 u ∂ x 2 = 0 .

pdepe требует по крайней мере одного параболического уравнения в системе. Другими словами, по крайней мере одно уравнение в системе должно включать производную времени.

pdepe также решает определенные 2D и 3-D задачи, которые уменьшают до 1D проблем из-за угловой симметрии (см. описание аргумента для симметрии постоянный m для получения дополнительной информации.

Partial Differential Equation Toolbox™ расширяет эту функциональность к обобщенным проблемам в 2D и 3-D с Дирихле и Неймановыми граничными условиями.

Решение 1D УЧП

1D УЧП включает функциональный u (x, t) , который зависит от времени t и одна пространственная переменная x. Решатель УЧП MATLAB pdepe решает системы 1D параболических и эллиптических УЧП формы

c ( x , t , u , ∂ u ∂ x ) ∂ u ∂ t = x − m ∂ ∂ x ( x m f ( x , t , u , ∂ u ∂ x ) ) + s ( x , t , u , ∂ u ∂ x ) .

Уравнение имеет свойства:

  • УЧП содержат для t0ttf и axb .
  • Пространственный интервал [a, b] должен быть конечным.
  • m может быть 0, 1, или 2, соответствуя плите , цилиндрической , или сферической симметрии, соответственно. Если m> 0 , то a ≥ 0 должен также содержать.
  • Коэффициент f ( x , t , u , ∂ u ∂ x ) термин потока и s ( x , t , u , ∂ u ∂ x ) характеристики выброса.
  • Термин потока должен зависеть от частной производной ∂u / ∂ x .

Связь частных производных относительно времени ограничивается умножением диагональной матрицей c ( x , t , u , ∂ u ∂ x ) . Диагональными элементами этой матрицы является или нуль или положительный. Элемент, который является нулем, соответствует эллиптическому уравнению, и любой другой элемент соответствует параболическому уравнению. Должно быть по крайней мере одно параболическое уравнение. Элемент c, который соответствует параболическому уравнению, может исчезнуть в изолированных значениях x, если они — точки mesh (точки, где решение оценено). Разрывы в c и s из-за материальных интерфейсов разрешены при условии, что точка mesh помещается в каждый интерфейс.

Процесс решения

Решить УЧП с pdepe , необходимо задать коэффициенты уравнения для c, f, и s, начальных условий, поведения решения на контурах и сетки точек, чтобы оценить решение на. Вызов функции sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) использование эта информация, чтобы вычислить решение на заданную mesh:

  • m постоянная симметрия.
  • pdefun определяет решаемые уравнения.
  • icfun задает начальные условия.
  • bcfun задает граничные условия.
  • xmesh вектор из пространственных значений для x.
  • tspan вектор из временных стоимостей для t.

Вместе, xmesh и tspan векторы формируют 2D сетку это pdepe оценивает решение на.

Уравнения

Необходимо описать УЧП в стандартной форме, ожидаемой pdepe . Написанный в этой форме, можно прочитать значения коэффициентов c, f и s.

В MATLAB можно закодировать уравнения с функцией формы

function [c,f,s] = pdefun(x,t,u,dudx) c = 1; f = dudx; s = 0; end

В этом случае pdefun определяет уравнение ∂ u ∂ t = ∂ 2 u ∂ x 2 . Если существует несколько уравнений, то c, f и s являются векторами с каждым элементом, соответствующим одному уравнению.

Начальные условия

В начальное время t = t 0, для всего x, компоненты решения удовлетворяют начальным условиям формы

u ( x , t 0 ) = u 0 ( x ) .

В MATLAB можно закодировать начальные условия с функцией формы

function u0 = icfun(x) u0 = 1; end

В этом случае u0 = 1 задает начальное условие u 0 (x, t 0) = 1 . Если существует несколько уравнений, то u0 вектор с каждым элементом, задающим начальное условие одного уравнения.

Граничные условия

За пределами x = a или x = b, для всего t, компоненты решения удовлетворяют граничным условиям формы

p ( x , t , u ) + q ( x , t ) f ( x , t , u , ∂ u ∂ x ) = 0.

q (x, t) является диагональной матрицей с элементами, которые являются или нулем или никогда не обнуляют. Обратите внимание на то, что граничные условия описываются в терминах потока f, а не частная производная u относительно x. Кроме того, этих двух коэффициентов p (x, t, u) и q (x, t) , только p может зависеть от u.

В MATLAB можно закодировать граничные условия с функцией формы

function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t) pL = uL; qL = 0; pR = uR - 1; qR = 0; end 

pL и qL коэффициенты для левого контура, в то время как pR и qR коэффициенты для правильного контура. В этом случае bcfun задает граничные условия

u L ( x L , t ) = 0 u R ( x R , t ) = 1

Если существует несколько уравнений, то выходные параметры pL , qL , pR , и qR векторы с каждым элементом, задающим граничное условие одного уравнения.

Опции интегрирования

Свойства интегрирования по умолчанию в решателе УЧП MATLAB выбраны, чтобы решить типичные проблемы. В некоторых случаях можно улучшать производительность решателя путем переопределения этих значений по умолчанию. Для этого использовать odeset создать options структура. Затем передайте структуру pdepe как последний входной параметр:

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

Из опций для базового решателя ОДУ ode15s , только показанные в следующей таблице доступны для pdepe .

Оценка решения

После того, как вы решаете уравнение с pdepe , MATLAB возвращает решение как трехмерный массив sol , где sol(i,j,k) содержит k компонент th решения оценен в t(i) и x(j) . В общем случае можно извлечь k компонент решения th с командой u = sol(. k) .

Mesh времени, которую вы задаете, используется просто для выходных целей и не влияет на внутренние временные шаги, взятые решателем. Однако пространственная mesh, которую вы задаете, может влиять на качество и скорость решения. После решения уравнения можно использовать pdeval оценивать структуру решения, возвращенную pdepe с различной пространственной mesh.

Пример: уравнение тепла

Примером параболического УЧП является уравнение тепла в одной размерности:

∂ u ∂ t = ∂ 2 u ∂ x 2 .

Это уравнение описывает рассеяние тепла для 0 ≤ x ≤ L и t ≥ 0 . Цель состоит в том, чтобы решить для температуры u ( x , t ) . Температура является первоначально ненулевой константой, таким образом, начальное условие

Кроме того, температура является нулем на левом контуре, и ненулевой на правильном контуре, таким образом, граничные условия

Чтобы решить это уравнение в MATLAB, необходимо закодировать уравнение, начальные условия и граничные условия, затем выбрать подходящую mesh решения прежде, чем вызвать решатель pdepe . Или можно включать необходимые функции как локальные функции в конце файла, (как в этом примере) или сохранить их как отдельные, именованные файлы в директории на пути MATLAB.

Уравнение кода

Прежде чем можно будет закодировать уравнение, необходимо убедиться что именно в форме pdepe решатель ожидает:

c ( x , t , u , ∂ u ∂ x ) ∂ u ∂ t = x — m ∂ ∂ x ( x m f ( x , t , u , ∂ u ∂ x ) ) + s ( x , t , u , ∂ u ∂ x ) .

В этой форме уравнение тепла

1 ⋅ ∂ u ∂ t = x 0 ∂ ∂ x ( x 0 ∂ u ∂ x ) + 0 .

Таким образом, значения коэффициентов следующие:

  • m = 0
  • c = 1
  • f = ∂ u ∂ x
  • s = 0

Значение m передается в качестве аргумента pdepe , в то время как другие коэффициенты закодированы в функции для уравнения, которое является

function [c,f,s] = heatpde(x,t,u,dudx) c = 1; f = dudx; s = 0; end 

(Примечание: Все функции включены как локальные функции в конце примера.)

Условие начальной буквы кода

Начальная функция условия для уравнения тепла присваивает постоянное значение для u 0 . Эта функция должна принять вход для x , даже если это не использовано.

function u0 = heatic(x) u0 = 0.5; end 

Граничные условия кода

Стандартная форма для граничных условий ожидается pdepe решатель

p ( x , t , u ) + q ( x , t ) f ( x , t , u , ∂ u ∂ x ) = 0 .

Написанный в этой форме, граничные условия для этой проблемы

u ( 0 , t ) + ( 0 ⋅ f ) = 0 ,

( u ( L , t ) — 1 ) + ( 0 ⋅ f ) = 0 .

Так значения для p и q

  • p L = u L , q L = 0 .
  • p R = u R — 1 , q R = 0 .

Соответствующая функция затем

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur - 1; qr = 0; end 

Выберите Solution Mesh

Используйте пространственную сетку 20 точек и сетку времени 30 точек. Поскольку решение быстро достигает устойчивого состояния, моменты времени рядом t = 0 более близко расположены вместе, чтобы получить это поведение в выходе.

L = 1; x = linspace(0,L,20); t = [linspace(0,0.05,20), linspace(0.5,5,10)];

Решите уравнение

Наконец, решите уравнение с помощью симметрии m , уравнение PDE, начальное условие, граничные условия и сетки для x и t .

m = 0; sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);

Постройте решение

Используйте imagesc визуализировать матрицу решения.

colormap hot imagesc(x,t,sol) colorbar xlabel('Distance x','interpreter','latex') ylabel('Time t','interpreter','latex') title('Heat Equation for $0 \le x \le 1$ and $0 \le t \le 5$','interpreter','latex')

Figure contains an axes object. The axes object with title Heat Equation for 0 less equals x less equals 1 and 0 less equals t less equals 5 contains an object of type image.

Локальные функции

function [c,f,s] = heatpde(x,t,u,dudx) c = 1; f = dudx; s = 0; end function u0 = heatic(x) u0 = 0.5; end function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur - 1; qr = 0; end

Примеры УЧП и файлы

Несколько доступных файлов в качестве примера служат превосходными начальными точками для наиболее распространенных 1D проблем УЧП. Чтобы исследовать и запустить примеры, используйте приложение Дифференциальных уравнений В качестве примера. Чтобы запустить это приложение, ввести

odeexamples

Чтобы открыть отдельный файл для редактирования, ввести

edit exampleFileName.m

Чтобы запустить пример, ввести

exampleFileName

Эта таблица содержит список доступных файлов УЧП в качестве примера.

Файл в качестве примера

Простой УЧП, который иллюстрирует формулировку, расчет и графический вывод решения.

Проблема, которая включает разрывы.

Проблема, которая требует вычисления значений частной производной.

Система двух УЧП, решение которых имеет пограничные слои и в концах интервала и изменяется быстро для маленького t.

Система УЧП со ступенчатыми функциями как начальные условия.

Ссылки

[1] Skeel, R. D. и М. Берзиньш, «Метод для Пространственной Дискретизации Параболических уравнений в Одной Пространственной переменной», SIAM Journal на Научном и Статистическом Вычислении, Издании 11, 1990, стр 1–32.

Смотрите также

Похожие темы

  • Решите один УЧП
  • Решите систему УЧП

Открытый пример

У вас есть модифицированная версия этого примера. Вы хотите открыть этот пример со своими редактированиями?

Как решать дифференциальные уравнения в матлабе

Теоретическая справка

Уравнения вида

решается путем n – кратного интегрирования.

не содержащее искомой функции y, подстановкой

где- низшая из производных, сводится к уравнению

порядок которого равен .

не содержащее независимой переменной , также допускает понижение порядка с помощью подстановки

Решение (аналитическое)

В левой части этого уравнения стоит третья производная искомой функции, в правой – функция только от , это уравнение вида (1.1).

Т.к. , то , , откуда

Поскольку , то последнее уравнение можно переписать так:

Интегрируя еще раз, находим общее решение исходного уравнения

Решение в среде MATLAB.

(Синтаксис команд MATLAB будет приведен ниже.)

1.Решим уравнение аналитически

Как видно, ответ совпадает с точностью до названия переменной.

2.Решим уравнение численно и построим его график.

Для этого необходимо составить задачу Коши (т.е. задать начальные условия),

а также интервал интегрирования и построения графика.

  • Создаем m.- файл, где задаем условие уравнения:
  • Задаем интервал интегрирования, начальные условия, и производим поиск решения.

Решим уравнение аналитически

Это уравнение вида (1.2) , для которого , .

Принимая низшую производную за новую неизвестную функцию , осуществим подстановку (1.3):

откуда . Т.о., данное уравнение сводится к уравнению первого порядка относительно :

из которого получаем

в котором правая часть зависит только от , т.е. уравнение вида(1.1).

Трижды интегрируя, получаем соответственно:

Решение в среде MATLAB

1.Решим уравнение аналитически

2.Решим уравнение численно и построим его график.

Для этого необходимо составить задачу Коши (т.е. задать начальные условия),

а также интервал интегрирования и построения графика.

Создаем m.- файл, где задаем условие уравнения:

Для избежания деления на 0 и “зацикливания” программы, добавляем к t малую величину 0,0000001, не влияющую, в общем, на численное решение.

Задаем интервал интегрирования, начальные условия, и производим поиск решения.

Решим уравнение аналитически

Это уравнение вида (1.5). Положим , тогда , поэтому уравнение имеет вид

Разделяя переменные (в предположении, что ) и интегрируя, получим

Замечание. Если , т.е. , то , , . Это решение получается из формулы (А) при .

Решение в среде MATLAB

1.Решим уравнение аналитически

2.Решим уравнение численно и построим его график.

Для этого необходимо составить задачу Коши (т.е. задать начальные условия),

а также интервал интегрирования и построения графика.

Создаем m.- файл, где задаем условие уравнения:

2 Линейные однородные уравнения с постоянными коэффициентами

Теоретическая справка

где — постоянные величины, называется линейным однородным уравнением — го порядка с постоянными коэффициентами.

где — его линейно независимые частные решения. Последние находятся с помощью характеристического уравнения

Если характеристическое уравнение имеет действительных различных корней , то каждому из них соответствует частное решение

и общее решение уравнения (2.1) принимает вид

Если уравнение (2.3) имеет действительных равных корней (т.е. корень имеет кратность ), то в формуле (2.2) им соответствуют решения

Однократным комплексно сопряженным корням уравнения (2.3) в формуле (2.2) соответствуют решения:

Комплексно сопряженным корням кратности соответствуют частные решения:

Решение (аналитическое)

Для данного уравнения с постоянными коэффициентами составим характеристическое уравнение( при этом нужно сохранить коэффициенты, вместо поставить 1, вместо ее производной -того порядка поставить):

Преобразуя правую часть уравнения , получим

В соответствии с формулой (2.5) получаем общее решение

Решение в MATLAB

Моделирование в MATLAB/Simulink и SCILAB/Scicos — page 231

дополнительных условий в дифференциальных уравнениях различают:

– все дополнительные условия заданы в одной (чаще начальной) точке

краевую задачу

– дополнительные условия указаны на границах

Большое количество уравнений может быть решено точно. Однако есть

уравнения, а особенно системы уравнений, для которых точное решение

записать нельзя. Такие уравнения и системы решают при помощи численных

методов. Так же численные методы применяют, если для уравнений с

известным аналитическим решением требуется найти числовое значение при

определенных исходных данных.

25. 2. Решатели ОДУ

Для решения дифференциальных уравнений и систем в Scilab

предусмотрена функция:

SEO Version

Warning.

You are currently viewing the SEO version of !text.
It has a number of design and functionality limitations.

We recommend viewing the Flash version or the basic HTML version of this publication.

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *