Глушко А.В., Глушко В.П. Современное программное обеспечение в пользовательском процессе: Сборник заданий по курсу

Подождите немного. Документ загружается.

Математический факультет

афедра уравнений в частных производных и теории вероятностей

Сборник за

аний

по курсу

" Современное программное

\n \n обеспечение в образовательном процессе "

Разделы :"Приближённое

решение дифференциальных уравнений",

"Численное решений

дифференциальных уравнений"

Null

для студентов второго курса дневного

отделения математического факультета

Составители :

А.В.Глушко,

В.П.Глушко

ВОРОНЕЖ 2002

ad ok7bis.nb 1

Предисловие

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

задач для обыкновенных дифференциальных уравнений (и систем

уравнений ) на ПК в среде Mathematica 4.1. Предлагаемые задания были

опробованы в 2001/ 2002 учебном году на 2-4 курсах математического

факультета ВГУ.

Каждое задание содержит пять задач. В первых двух задачах

предлагается решить задачу Коши для дифференциального уравнения

первого порядка, причём во второй задаче требуется обеспечить

повышенную точность (по сравнению со стандартной машинной

точностью). В третьей задаче рассматривается система из трёх

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

график решения этой системы в трёхмерном фазовом пространстве.

В четвёртой задаче изучается серия из пяти начальных задач для

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

Предлагается построить не только графики решений каждой из задач, но

и объединённый график решений, позволяющий провести

сравнительный анализ изменений решений в зависимости от изменений

параметра.

В последней пятой задаче рассматривается граничная задача на

конечном отрезке вещественной оси для одного уравнения второго

порядка с комплекснозначными данными. Поскольку решение

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

Предлагается два способа графического изображения решения в этом

случае.

Отметим также, что программа Mathematica 4.1 обладает несколькими

методами численного решения дифференциальных уравнений. По

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

ассмотренных задач авторы предпочли метод Рунге-Кутта.

Решение каждой задачи сопровождается его проверкой и выяснением

фактической точности построенного приближённого решения. Для этого

предлагается простая оригинальная методика.

В первой части сборника приведён пример выполнения задания с

соответствующими пояснениями. Более полные комментарии к

ad ok7bis.nb 2

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

указанных в конце сборника в списке рекомендуемой литературы, а

также непосредственно в справочнике HelpBrowser (на английском

языке) в самом пакете Mathematica 4.1.

Пример выполнения задания

Задача 1

Найти численное решение следующей начальной задачи для одного

уравнения первого порядка. Построить график решения и провести

проверку

8- Cos@ x

@xDD +

@xD == Sin@x

@0D ä 1<, 8x, -5, 5<

Введём данные задачи

f10 = Sin@x

D; w10 = 1;

eq10 = 8- Cos@x * #@xDD + D@#@xD,xD,#@xD

.xÆ 0< &

8-Cos@x#1@xDD+∑

#1@xD,#1@xD

.xØ 0< &

Теперь задачу можно записать в виде

g10 = 88eq10@#D@@1DD == f10, eq10@#D@@2DD == w10<, 8x, -5, 5<< &

88eq10@#1DP1T == f10, eq10@#1DP2T == w10<, 8x, -5, 5<< &

Для проверки набираем

g10@yD

88-Cos@xy@xDD+ y

@xD == Sin@x

D,y@0D == 1<, 8x, -5, 5<<

После этого находим решение задачи

ad ok7bis.nb 3

y10 = y ê. First@NDSolve@g10@yD@@1DD, y, g10@yD@@2DDDD

InterpolatingFunction@88-5., 5.<<, <>D

pl10 = Plot@y10@xD, 8x, -5, 5<, PlotPoints Æ 100, AxesLabel Æ 8"x", "y"<D

-4 -2 2 4

-1

-0.5

0.5

1.5

Ü Graphics Ü

Перейдём к проверке построенного решения. Вначале проверим выполнение

начального условия

y10@0D - w10

Для проверки выполнения уравнения будем подсчитывать разность между

левой и правой частями уравнения

eq10@yD@@1DD == f10

в точках (в данном случае 500), выбранных случайным образом на

интервале (-5,5). Затем после округления (с помощью команды Chop)

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

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

может занять довольно много времени (которое зависит, конечно, от типа

ПК), следует для контроля включить счётчик времени выполнения команды.

<< Utilities`ShowTime`

ad ok7bis.nb 4

Max@Abs@Table@x10 = Random@Real, 8-5, 5<D;

Chop@HHeq10@y10D@@1DD - f10L

.xÆ x10L,1* 10^-4D, 8500<DDD

57.34 Second

0.000197888

Счётчик времени после использования желательно выключить.

Off@ShowTimeD

0. Second

Задача 2

Найти численное решение с повышенной точностью следующей

начальной задачи для одного уравнения первого порядка. Построить

график решения и провести проверку

8Tanh@ x

* y@xD

D + y

@xD == Cos@x

D,y@0D == 0.5<, 8x, -3.5, 3.5<,

MaxSteps Æ 16000, AccuracyGoal -> •, PrecisionGoal Æ 18,

WorkingPrecision Æ 18, Method -> RungeKutta

Вначале введём данные задачи

f20 = Cos@x^3D; w20 = 0.5;

eq20 = 8D@#@xD,xD + Tanh@x^4* #@xD^2D,#@xD

.xÆ 0< &

8∑

#1@xD+ Tanh@x

#1@xD

D,#1@xD

.xØ 0< &

Запишем задачу в виде

g20 = 88eq20@#D@@1DD == f20, eq20@#D@@2DD == w20<, 8x, -3.5, 3.5<< &

88eq20@#1DP1T == f20, eq20@#1DP2T == w20<, 8x, -3.5, 3.5<< &

Программа Mathematica по умолчанию (default) проводит все внутренние

вычисления с машинной точностью WorkingPrecisionØ16. В

подавляющем большинстве задач этой точности бывает достаточно. Однако

в ряде задач требуется повысить точность вычислений. Найдём решение

задачи g20, установив повышенную точность вычислений с помощью

опций: AccuracyGoal

→∞

, PrecisionGoal

→

18, Working-

Precision

→

ad ok7bis.nb 5

y20 = y ê. First@NDSolve@g20@yD@@1DD, y, g20@yD@@2DD,

MaxSteps Æ 16000 , AccuracyGoal Æ•, PrecisionGoal Æ 18,

WorkingPrecision Æ 18, Method Æ RungeKuttaDD

InterpolatingFunction@88-3.5, 3.5<<, <>D

pl20 = Plot@y20@xD, 8x, -3.5, 3.5<, PlotRange Æ 8-0.45, 1.25<,

PlotPoints Æ 120, AxesLabel Æ 8"x", "y"<D

-3 -2 -1 1 2 3

-0.25

0.25

0.5

0.75

1.25

Ü Graphics Ü

Если часть графика функции выглядит слишком мелкой, то её можно

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

pl201 = Plot@y20@xD, 8x, -3.5, -2.2<, PlotRange Æ 8-0.1, 0.2<,

PlotPoints Æ 120, AxesLabel Æ 8"x", "y"<D

-3.4 -3.2 -2.8 -2.6 -2.4 -2.2

-0.1

-0.05

0.05

0.1

0.15

0.2

Ü Graphics Ü

Проведём проверку решения

ad ok7bis.nb 6

y20@0D - w20

<< Utilities`ShowTime`

Max@Abs@Table@x20 = Random@Real, 8-3.5, 3.5<D;

Chop@HHeq20@y20D@@1DD - f20L

.xÆ x20L,1* 10^-6D, 8500<DDD

115.13 Second

2.86454 × 10

−6

Счётчик времени после использования следует выключить.

Off@ShowTimeD

Задача 3

Найти численное решение следующей начальной задач для

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

график решения в фазовом пространстве с помощью команды Paramet-

ricPlot3D в подпакете <<Graphics`ParametricPlot3D`. Провести проверку

решения.

8-H2 + 0.05 Sin@tDL Hy@tD + z@tDL + x

@tD == 0,

H1 - 0.06 Cos@tDL Hx@tD + z@tDL + y

@tD == 0,

H1 + 0.001 t

L x@tD - H1 + 0.07 tL y@tD

+ z

@tD == 0,

x@0D == 0.99, y@0D == 0.02, z@0D == -1.97<, 8t, 0, 4.15<,

MaxSteps Æ 15000, AccuracyGoal -> •, PrecisionGoal -> 17,

WorkingPrecision -> 17, Method -> RungeKutta

Начальные данные задачи имеют вид

w301 = 0.99; w302 = 0.02; w303 =-1.97;

Введём оператор

ad ok7bis.nb 7

eq30 = 8D@#1@tD,tD - H2 + 0.05 * Sin@tDL* H#2@tD + #3@tDL,

D@#2@tD,tD + H1 - 0.06 * Cos@tDL * H#1@tD + #3@tDL,

D@#3@tD,tD + H1 + 0.001 * t^2L* #1@tD - H1 + 0.07 * tL* #2@tD

#1@tDê.t-> 0, #2@tDê.t-> 0, #3@tD

.t-> 0< &

8∑

#1@tD- H2 + 0.05 Sin@tDLH#2@tD+ #3@tDL,

∑

#2@tD+ H1 - 0.06 Cos@tDLH#1@tD+ #3@tDL,

∑

#3@tD+ H1 + 0.001 t

L#1@tD- H1 + 0.07 tL#2@tD

#1@tDê.tØ 0, #2@tDê.tØ 0, #3@tD

.tØ 0< &

Для проверки набираем

eq30@x, y, zD

8-H2 + 0.05 Sin@tDLHy@tD+ z@tDL+ x

@tD, H1 - 0.06 Cos@tDLHx@tD+ z@tDL+ y

@tD,

H1 + 0.001 t

Lx@tD - H1 + 0.07 tLy@tD

+ z

@tD,x@0D,y@0D,z@0D<

Обозначим

g30 = 88eq30@#1, #2, #3D@@1DD ä 0, eq30@#1, #2, #3D@@2DD ä 0,

eq30@#1, #2, #3D@@3DD ä 0, eq30@#1, #2, #3D@@4DD ä w301,

eq30@#1, #2, #3D@@5DD ä w302, eq30@#1, #2, #3D@@6DD == w303<, 8t, 0, 4.15<< &

88eq30@#1, #2, #3DP1T == 0, eq30@#1, #2, #3DP2T == 0, eq30@#1, #2, #3DP3T == 0,

eq30@#1, #2, #3DP4T == w301, eq30@#1, #2, #3DP5T == w302,

eq30@#1, #2, #3DP6T == w303<, 8t, 0, 4.15<< &

Для проверки набираем

g30@x, y, zD

88-H2 + 0.05 Sin@tDLHy@tD+ z@tDL+ x

@tD == 0,

H1 - 0.06 Cos@tDLHx@tD+ z@tDL+ y

@tD == 0,

H1 + 0.001 t

Lx@tD - H1 + 0.07 tLy@tD

+ z

@tD == 0,

x@0D == 0.99, y@0D == 0.02, z@0D == -1.97<, 8t, 0, 4.15<<

Находим решение задачи

ad ok7bis.nb 8

8x30, y30, z30< =

8x, y, z<ê. First@NDSolve@g30@x, y, zD@@1DD, 8x, y, z<, g30@x, y, zD@@2DD,

MaxSteps Æ 15000, AccuracyGoal -> •, PrecisionGoal -> 17,

WorkingPrecision -> 17, Method -> RungeKuttaDD

8InterpolatingFunction@880, 4.15<<, <>D,

InterpolatingFunction@880, 4.15<<, <>D, InterpolatingFunction@880, 4.15<<, <>D<

Вызовем подпакет приложения <<Graphics`ParametricPlot3D`

<< Graphics`ParametricPlot3D`

и построим график решения в фазовом пространстве {x,y,z}

sty301 = 8Thickness@0.0180 - 0.002 * #D, GrayLevel@0.0 + 0.035 * #D< &;

sty302 = 8Hue@# ê 4.5D, Thickness@0.0150 − 0.00085 ∗ #D< &;

ad ok7bis.nb 9

ParametricPlot3D@Evaluate@8x30@tD, y30@tD, z30@tD, sty301@tD<D,

Evaluate@g30@x, y, zD@@2DDD, AspectRatio Æ Automatic,

AxesLabel Æ 8"x", "y", "z"<, BoxRatios Æ 82.5, 4.0, 3.0<,

ViewPoint -> 80.880, -2.380, 0.620<D

-1

-2

-1

 Graphics3D 

Этот же график в цвете строится по команде

ParametricPlot3D@Evaluate@8x30@tD, y30@tD, z30@tD, sty302@tD<D,

Evaluate@g30@x, y, zD@@2DDD, AspectRatio Æ Automatic,

AxesLabel Æ 8"x", "y", "z"<, BoxRatios Æ 82.5, 4.0, 3.0<,

ViewPoint -> 80.880, -2.380, 0.620<D

Проверка решения. Вначале проверяются начальные условия.

Chop@8eq30@x30, y30, z30D@@4DD - w301,

eq30@x301, y301, z301D@@5DD - w302,

eq30@x30, y30, z30D@@6DD - w303<,10^-15D

80, 0, 0<

Затем проверяется выполнение трёх уравнений

ad ok7bis.nb 10