И. И. Мечникова лаборатория кафедра компьютерных методов экспериментальной экспериментальной физики физики компьютерный практикум

Вид материала

Содержание

Работа №5 Решение систем обыкновенных дифференциальных уравнений 1. Краткие теоретические сведения
2. Реализация метода Рунге-Кутта 4-го порядка в библиотечном модуле Turbo-Pascal
3. Примеры использования процедуры Runge
В а р и а н т 1
В а р и а н т 2
В а р и а н т 3
В а р и а н т 4
В а р и а н т 5
В а р и а н т 6

Подобный материал:

ОДЕССКИЙ ГОСУНИВЕРСИТЕТ им. И.И.Мечникова

ЛАБОРАТОРИЯ	КАФЕДРА
КОМПЬЮТЕРНЫХ МЕТОДОВ	ЭКСПЕРИМЕНТАЛЬНОЙ
ЭКСПЕРИМЕНТАЛЬНОЙ ФИЗИКИ	ФИЗИКИ

КОМПЬЮТЕРНЫЙ ПРАКТИКУМ

ДЛЯ ФИЗИКОВ

МЕТОДИЧЕСКИЕ УКАЗАНИЯ ДЛЯ СТУДЕНТОВ

ФИЗИЧЕСКОГО ФАКУЛЬТЕТА

Автор доцент П.А.Виктор

Р а б о т а 5

РЕШЕНИЕ СИСТЕМ ОБЫКНОВЕННЫХ

ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

ОДЕССА - 1995

Работа №5

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

1. Краткие теоретические сведения

При математическом моделировании различных процессов, протекающих в тех или иных физических системах, часто приходится иметь дело с обыкновенными дифференциальными уравнениями или системами этих уравнений, в которых в качестве независимой переменной выступает время. Нахождение функций, являющихся решениями дифференциальных уравнений, аналитическими методами часто оказывается невозможным. В то же время, разработано множество численных методов, позволяющих находить приближенные решения практически с любой необходимой точностью. Одним из наиболее распространенных является метод Рунге-Кутта 4-го порядка, реализованный в описанной ниже стандартной процедуре Runge.

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

µ§ (1),

содержащее независимую переменную x, неизвестную функцию µ§ и ее производные µ§. Порядок уравнения определяется наивысшим порядком содержащихся в нем производных.

Функция µ§ называется решением дифференциального уравнения (1), если при подстановке ее в уравнение последнее обращается в тождество. Каждое уравнение имеет бесконечное множество решений. Для получения конкретного частного решения необходимо указать начальные условия, т.е. задать числовые значения функции и ее µ§ производных в некоторой фиксированной точке µ§:

µ§ (2).

Уравнение (1), записанное в виде

µ§ (3),

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

Совокупность дифференциальных уравнений, содержащих несколько неизвестных функций и их производных, составляет систему дифференциальных уравнений. В частности, система уравнений 1-го порядка, разрешенных относительно производных, имеет вид:

µ§ (4).

Здесь µ§ неизвестные функции общего аргумента x. Для нахождения частного решения системы (4) начальные условия задаются в виде

µ§ .

Известно, что любое уравнение вида (3) n-го порядка можно свести к эквивалентной системе n уравнений 1-го порядка вида (4) путем введения вспомогательных переменных µ§. Аналогичным образом к системе уравнений 1-го порядка можно свести систему нескольких уравнений более высоких порядков. Поэтому математическое обеспечение компьютеров предусматривает не решение одного уравнения высокого порядка, а решение системы нескольких уравнений 1-го порядка.

Например, уравнение 2-го порядка

µ§ (5)

с начальными условиями µ§ после разрешения относительно старшей производной и введения вспомогательных переменных µ§ сводится к системе

µ§ (6)

с начальными условиями µ§ . Здесь неизвестная в уравнении (5) функция µ§ есть решение µ§ системы (6), а решение µ§ представляет собой производную этой функции.

2. Реализация метода Рунге-Кутта 4-го порядка в библиотечном модуле Turbo-Pascal

Процедура интегрирования системы обыкновенных дифференциальных уравнений входит в библиотечный модуль Numerics под именем Runge. Ниже приводится полный текст той части модуля, которая относится к данной процедуре. Он может представлять интерес для составителей собственных библиотек численных методов и как пример использования в Pascal-программе таких современных средств, как процедурные типы, нетипизированне переменные и совмещение переменных в памяти. Пользователи могут пропустить материал до конца данного раздела.

Процедура вычисляет значения функций µ§ в точке µ§ по значениям µ§ в исходной точке x, по следующей схеме:

µ§

Здесь через µ§ обозначены совокупность искомых функций µ§ и функции правых частей системы µ§ , соответственно.

unit Numerics;

interface

{$N+,E+,F+}

const

ArraySize = (2*MaxInt) div SizeOf (extended);

type

ArrayType = array [1..ArraySize] of extended;

ProcType = procedure (T: extended; var YY, FF);

procedure Runge

{ pешение систем обыкновенных диффеpенциальных

уpавнений методом Рунге-Кутта 4-го поpядка }

(F : ProcType; { пpоцедуpа вычисления пpавых частей }

NEqn : integer; { количество диффеpенциальных уpавнений }

var Yi; { пpи входе - вектоp начальных условий,

пpи выходе - вектоp pешения }

var T:extended; { пpи входе - начальная точка

интегpиpования, пpи выходе T = Tout }

Tout : extended; { конечная точка интегpиpования }

h : extended; { шаг интегpиpования }

var Work); { pабочий массив длиной 3*NEqn }

{ ============================================== }

implementation

const

Eps: extended = 0; { нач. значение машинного эпсилон }

{ ---------------------------------------------- }

procedure Runge (F: ProcType; NEqn: integer; var Yi;

var T: extended;

Tout: extended; h: extended; var Work);

var

W: array [1..3, 1..ArraySize div 3]

of extended absolute Work;

Y: ArrayType absolute Yi;

k: integer;

dT: extended;

begin

{ пpовеpка пpавильности паpаметpов вызова }

if (NEqn <= 0) or (h = 0) then

begin

Writeln ('Невеpный паpаметp вызова Runge');

Halt

end;

{ вычисление машинного эпсилон пpи пеpвом вызове }

if Eps = 0 then

begin

Eps:= 1;

while 1+Eps/2 <> 1 do

Eps:= Eps/2;

Eps:= 20*Eps

end;

{ начало основного цикла }

repeat

{ опpеделяем шаг, чтобы точно попасть в Tout }

F (T, Y, W);

dT:= Tout - T;

if Abs (dT) <= Abs (T) * Eps then

{ T в Eps-окpестности Tout, экстpаполяция и выход }

begin

for k:= 1 to NEqn do

Y[k]:= Y[k] + dT*W[1,k];

T:= Tout;

Exit

end;

if Abs (dT) < Abs (h) then

h:= dT;

if dT < 0 then

h:= -Abs (h)

else

h:= Abs (h);

{ основные вычисления по фоpмулам Рунге-Кутта }

for k:= 1 to NEqn do

W[2,k]:= Y[k] + h*W[1,k]/2;

F (T+h/2, W[2,1], W[3,1]);

for k:= 1 to NEqn do

begin

W[1,k]:= W[1,k] + 2*W[3,k];

W[2,k]:= Y[k] + h*W[3,k]/2

end;

F (T+h/2, W[2,1], W[3,1]);

for k:= 1 to NEqn do

begin

W[1,k]:= W[1,k] + 2*W[3,k];

W[2,k]:= Y[k] + h*W[3,k]

end;

F (T+h, W[2,1], W[3,1]);

for k:= 1 to NEqn do

Y[k]:= Y[k] + h*(W[1,k]+W[3,k])/6;

T:= T + h

until T = Tout

end; { Runge }

END.

3. Примеры использования процедуры Runge

Процедура Runge для своей работы использует 6 параметров:

Runge (F, NEqn, Yi, T, Tout, h, Work);

Здесь

F - имя процедуры, в которой производится вычисление функций правых частей системы уравнений. Эта процедура должна иметь вполне определенную структуру (ей в модуле Numerics соответствует тип ProcType), которая будет описана ниже;

NEqn - число уравнений в системе (integer);

Yi - имя массива зависимых переменных типа extended, который перед вызовом процедуры Runge должен содержать начальные значения этих переменных, а после выполнения процедуры в этот массив будут записаны результаты решения системы;

T - независимая переменная (extended), которой перед вызовом процедуры должно быть присвоено ее начальное значение, а после выполнения T становится равным Tout;

Tout - конец интервала интегрирования (extended);

h - шаг интегрирования (extended);

Work - вспомогательный массив переменных типа extended размером 3*NEqn, в котором процедура Runge размещает промежуточные результаты вычислений.

Процедура вычисления правых частей уравнений должна быть построена по следующему шаблону:

procedure Right (X: extended; var YY, FF);

var

Y: array [1..NEqn] of extended absolute YY;

F: array [1..NEqn] of extended absolute FF;

begin

F[1]:= .... {правая часть первого уравнения};

F[2]:= .... {правая часть второго уравнения};

....

end; { Right }

Здесь массив зависимых переменных Y совмещается в памяти с так называемой нетипизированной переменной YY с помощью зарезервированного слова absolute (т.е. одна и та же область памяти закрепляется как за массивом Y , так и за переменной YY). Поскольку YY является параметром процедуры, то такой механизм дает возможность эту же область памяти использовать и другим процедурам (в нашем случае Runge). Так из одной процедуры в другую могут передаваться массивы произвольной длины.

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

Следует помнить, что поскольку процедура Runge использует переменные типа extended, то в вызывающей ее программе необходимо предусмотреть подключение математического сопроцессора, а использование процедурных типов (ProcType) требует включения механизма дальнего вызова подпрограмм (см. работу №4). Это достигается помещением в начало программы директивы компилятору вида {$N+,E+,F+} .

Рассмотрим ряд примеров использования процедуры Runge.

А) Простейший случай.

Найдем решение системы (6) с указанными начальными условиями в точке x=0.4 с точностью 0.0001 . Данную задачу решает программа DIF1. Здесь же приведены результаты ее работы.

program DIF1;

{$N+,E+,F+}

uses Numerics;

const

N = 2; { число уpавнений в системе }

var

Y: array [1..N] of extended; { зависимые пеpеменные }

X: extended; { независимая пеpеменная }

Work: array [1..3*N] of extended; { pабочий массив }

procedure Right (X: extended; var YY, FF);

var

Y: array [1..N] of extended absolute YY;

F: array [1..N] of extended absolute FF;

begin

F[1]:= Y[2];

F[2]:= (2*Y[1] - Sqr(Y[2]) - 1) / (Sqr(X) + 1)

end; { Right }

BEGIN

{ задание начальных условий }

X:= 0; Y[1]:= 0; Y[2]:= 1;

{ вызов пpоцедуpы Runge }

Runge (Right, N, Y, X, 0.4, 0.1, Work);

{ вывод pезультатов }

Writeln (X:5:1, Y[1]:15:6, Y[2]:15:6);

Readln

END.

0.4 0.288091 0.530858

Программа начинается с директив подключения математического сопроцессора и механизма дальнего вызова подпрограмм. Затем следует объявление об использовании библиотечного модуля Numerics.

Число уравнений в системе удобнее всего задавать в виде константы, которая в дальнейшем используется в описаниях всех массивов. Это прежде всего массив Y с начальными значениями зависимых переменных, на место которых процедура Runge запишет решения системы, и рабочий массив Work. Кроме того, это массивы Y и F процедуры Right.

Тело процедуры Right составляют формулы вычисления правых частей решаемой системы уравнений.

Раздел операторов программы начинается с задания начальных условий, затем следует вызов процедуры Runge с шагом интегрирования 0.1 и вывод на экран значений независимой переменной X и решений системы Y[1] и Y[2].

Для оценки точности вычислений повторим расчет, заменив величину шага на 0.05 . При этом получим:

0.4 0.288098 0.530861

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

Б) Получение таблицы решений системы.

Как отмечалось выше, после выполнения процедуры Runge результаты решения системы записываются на место начальных условий. Следовательно, эти результаты можно использовать в качестве начальных условий для нахождения решений в следующей точке. Многократное использование такого приема в циклической программе позволяет получить таблицу решений системы в некотором интервале изменения аргумента с заданным шагом (шагом печати). Внутри цикла должно происходить увеличение верхнего предела интегрирования на величину шага печати.

Описанный алгоритм реализован в программе DIF2, построенной на основе программы DIF1. Решение системы (6) проводится в интервале изменения аргумента X от 0.2 до 1 с шагом печати 0.2 . Шаг интегрирования 0.1 .

program DIF2;

{$N+,E+,F+}

uses Numerics;

const

N = 2; { число уpавнений в системе }

var

Y: array [1..N] of extended; { зависимые пеpеменные }

X: extended; { независимая пеpеменная }

Work: array [1..3*N] of extended; { pабочий массив }

i: integer;

procedure Right (X: extended; var YY, FF);

{ вычисление пpавых частей системы }

var

Y: array [1..N] of extended absolute YY;

F: array [1..N] of extended absolute FF;

begin

F[1]:= Y[2];

F[2]:= (2*Y[1] - Sqr(Y[2]) - 1) / (Sqr(X) + 1)

end; { Right }

BEGIN

Writeln ('X':4, 'Y1':14, 'Y2':15); { заголовок таблицы }

Writeln;

{ задание начальных условий }

X:= 0; Y[1]:= 0; Y[2]:= 1;

{ цикл интегpиpования системы }

for i:= 1 to 5 do

begin

Runge (Right, N, Y, X, 0.2*i, 0.1, Work);

Writeln (X:5:1, Y[1]:15:6, Y[2]:15:6)

end;

Readln

END.

X Y1 Y2

0.2 0.166858 0.698095

0.4 0.288091 0.530858

0.6 0.384320 0.441231

0.8 0.467633 0.397784

1.0 0.545280 0.382282

Обратим внимание на то, что конец текущего интревала интегрирования выражается через через параметр цикла i непосредственно при обращении к процедуре Runge. Поэтому для него не требуется отдельного идентификатора.

В) Случай скачкообразного изменения внешних условий.

При моделировании физических процессов, протекающих в условиях, меняющихся скачкообразно, приходится иметь дело с дифференциальными уравнениями или их системами, правые части которых при разных значениях независимой переменной x описываются различными функциями. Например, пусть в системе уравнений

µ§ правые части имеют вид

µ§ .

Тогда функции µ§ удобно представить в форме

µ§ .

Программа решения системы должна предусматривать присваивание А значения 0 или 1 в зависимости от величины x. Для правильного изменения значения A точка µ§ должна являться границей одного из текущих промежутков интегрирования.

В качестве примера рассмотрим программу решения уравнения первого порядка

µ§ .

с начальными условиями µ§ (программа DIF3). Решение проводится в интервале µ§ с шагом печати результатов 0.5 .

program DIF3;

{$N+,E+,F+}

uses Numerics;

const

N = 1; { число уpавнений в системе }

var

Y: array [1..N] of extended; { зависимые пеpеменные }

X: extended; { независимая пеpеменная }

Work: array [1..3*N] of extended; { pабочий массив }

A,i: integer;

procedure Right (X: extended; var YY, FF);

{ вычисление пpавых частей системы }

var

Y: array [1..N] of extended absolute YY;

F: array [1..N] of extended absolute FF;

begin

F[1]:= 10*A - Y[1]

end; { Right }

BEGIN

Writeln ('X':4, 'Y1':14); { заголовок таблицы }

Writeln;

{ задание начальных условий }

X:= 0; Y[1]:= 0;

{ цикл интегpиpования системы }

A:= 1;

for i:= 1 to 4 do

begin

Runge (Right, N, Y, X, 0.5*i, 0.1, Work);

Writeln (X:5:1, Y[1]:15:6)

end;

A:= 0;

for i:= 5 to 8 do

begin

Runge (Right, N, Y, X, 0.5*i, 0.1, Work);

Writeln (X:5:1, Y[1]:15:6)

end;

Readln

END.

X Y1

0.5 3.934691

1.0 6.321202

1.5 7.768695

2.0 8.646645

2.5 5.244457

3.0 3.180926

3.5 1.929330

4.0 1.170198

Здесь правая часть уравнения преобразована к виду

µ§ .

Основу программы составляют два следующих один за другим цикла for..to..do . Перед первым циклом, где интегрирование проводится до точки "скачка" x=2, параметру А присваивается единичное значение. Перед вторым циклом выполняется присвавание этому параметру нулевого значения.

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

program DIF4;

{$N+,E+,F+}

uses Numerics;

var

Y: extended; { зависимая пеpеменная }

X: extended; { независимая пеpеменная }

Work: array [1..3] of extended; { pабочий массив }

A,i: integer;

OutFile: text;

procedure Right (X: extended; var YY, FF);

{ вычисление пpавой части уpавнения }

var

Y: extended absolute YY;

F: extended absolute FF;

begin

F:= 10*A - Y

end; { Right }

BEGIN

Assign (OutFile, 'Dif4.dat');

Rewrite (OutFile);

Writeln (OutFile, 'X':4, 'Y1':14); { заголовок таблицы }

Writeln (OutFile);

{ задание начальных условий }

X:= 0; Y:= 0;

{ цикл интегpиpования уpавнения }

for i:= 1 to 8 do

begin

if i <= 4 then

A:= 1

else

A:= 0;

Runge (Right, 1, Y, X, 0.5*i, 0.1, Work);

Writeln (OutFile, X:5:1, Y:15:6)

end;

Close (OutFile)

END.

Данная программа содержит только один цикл for..to..do , внутри которого изменение А с 1 на 0 происходит при выполнении условия µ§. Отметим также, что в программе вместо массивов Y[1..1] и F[1..1], состоящих из одного элемента (решается всего одно уравнение, N=1) используются простые переменные Y и F. Кроме того, результаты решения уравнения выводятся не на экран, а в текстовый файл DIF4.DAT .

________________________

4. Задание

1. В соответствии с одним из предложенных ниже вариантов напишите программу интегрирования системы дифференциальных уравнений методом Рунге-Кутта. Предусмотрите вывод результатов расчета в виде таблицы, записываемой в дисковый файл с расширением .DAT.

2. Проведите расчет по составленной программе. Пользуясь полученным файлом с результатами расчета, постройте графики исследуемых зависимостей либо вручную на миллиметровой бумаге, либо с помощью сервисной программы GRAPHER.

В а р и а н т 1

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

µ§ , где

j - угол отклонения от положения равновесия, w - угловая частота малых колебаний. Рассчитайте зависимости µ§, в интервале µ§ с шагом 0.5 с для µ§.

Начальные условия: µ§. Шаг интегрирования примите равным 0.1 .

В а р и а н т 2

Если в некотором сосуде в результате распада радиоактивного препарата А образуется радиоактивный препарат В, который, распадаясь, превращается в стабильный продукт С, то количества этих препаратов µ§ могут быть рассчитаны путем решения следующей системы обыкновенных дифференциальных уравнений 1-го порядка:

µ §где µ§ постоянные распада препаратов А и В, соответственно.

Рассчитайте зависимости µ§, µ§ и µ§ в интервале от 0 до 20 с с шагом 1 с, полагая

µ§

Начальные условия: µ§ Шаг интегрирования 0.5 .

В а р и а н т 3

Вынужденные колебания гармонического осциллятора описываются уравнением

µ§ , где

w - угловая частота, d - коэффициент затухания, m - масса осциллятора, µ§ - вынуждающая сила. Рассчитайте зависимость µ§ в интервале µ§ с шагом 0.075с , если колебания возбуждаются под действием силы

µ§ .

Начальные условия: µ§. Пусть µ§, µ§. Шаг интегрирования 0.005 .

В а р и а н т 4

µ §

Движение тела 3 в гравитационном поле неподвижных тел 1 и 2 (расположение тел показано на рисунке) описывается уравнениями

µ§

где µ§ гравитационная постоянная, µ§ массы тел 1 и 2, соответственно. Рассчитайте зависимости µ§ и µ§ в интервале µ§ с шагом 2, полагая a=4, b=10, d=10. Постройте траекторию движения тела 3 µ§.

Начальные условия: µ§.

Шаг интегрирования 0.5 .

В а р и а н т 5

Уравнение движения снаряда, выпущенного из орудия под углом к горизонту, с учетом сил сопротивления воздуха имеет вид:

µ§

µ§ проекции скорости на координатные оси,

µ§ модуль скорости,

µ§ ускорение свободного падения,

µ§ масса снаряда,

µ§ коэффициент аэродинамического сопротивления,

µ§ площадь поперечного сечения снаряда,

µ§ плотность воздуха.

Рассчитайте траектории полета снаряда µ§ для промежутка времени µ§ с шагом 1с

а) считая плотность воздуха постоянной µ§;

б) учитывая зависимость плотности воздуха от высоты

µ§ .

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

µ§ ,

шаг интегрирования 0.5 с .

В а р и а н т 6

Процесс изменения напряжений µ§ на конденсаторах µ§ в изображенном на рисунке фазовращателе описывается системой уравнений

µ § µ§

где µ§ напряжение на входе фазовращателя. Рассчитайте зависимости µ§, если в момент t=0 на фазовращатель было подано на время 4с постоянное напряжение µ§. Пусть µ§, µ§.

Вычисления проведите в интервале µ§ с шагом 0.25с . Начальные условия: µ§. Шаг интегрирования 0.125 .

* * *

Blog