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

КУРСОВАЯ РАБОТА

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

• ВВЕДЕНИЕ
• 1. ЧИСЛЕННОЕ РЕШЕНИЕ МЕТОДОМ ЭЙЛЕРА
• 1.1 Постановка задачи
• 1.2 Математическая модель задачи (метод Эйлера)
• 1.3 Исходные данные
• 1.4 Численное решение уравнения методом Эйлера в Excel
• 2. ЧИСЛЕННОЕ РЕШЕНИЕ МЕТОДОМ РУНГЕ-КУТТА В EXCEL И TURBO PASCAL 7.0
• 2.1 Постановка задачи
• 2.2 Математическая модель задачи (метод Рунге-Кутта)
• 2.3 Исходные данные
• 2.4 Численное решение уравнения методом Рунге-Кутта в Excel
• 2.5 Блок-схема алгоритма
• 2.6 Программа на языке Turbo Pascal
• 2.7 Выполнение расчетов
• 2.8 Результаты расчетов
• 2.9 Представление результатов в виде графиков
• 2.10 Анализ результатов
• 3. МЕТОД РУНГЕ-КУТТА ДЛЯ ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ ВТОРОГО ПОРЯДКА
• 3.1 Постановка задачи
• 3.2 Математическая модель задачи
• 3.3 Исходные данные
• 3.4 Расчет контрольного варианта в Excel для тестирования
• 3.5 Блок-схема алгоритма
• 3.6 Численное решение задачи с использованием Turbo Pascal 7.0.
• 3.7 Выполнение расчетов
• 3.8 Результаты расчетов
• 3.9 Представление результатов в виде графиков
• 3.10 Анализ результатов
• 4. МОДЕЛЬ ТИПА «ХИЩНИК-ЖЕРТВА» С УЧЕТОМ ВНУТРИВИДОВОГО ВЗАИМОДЕЙСТВИЯ
• 4.1 Постановка задачи
• 4.2 Математическая модель задачи
• 4.3 Исходные данные
• 4.4 Расчет контрольного варианта в Excel для тестирования
• 4.5 Блок-схема алгоритма
• 4.6 Численное решение задачи с использованием Turbo Pascal 7.0.
• 4.7 Выполнение расчетов
• 4.8 Результаты расчетов
• 4.9 Представление результатов в виде графиков
• 4.10 Анализ результатов
• ЗАКЛЮЧЕНИЕ
• СПИСОК ИСПОЛЬЗОВАННОЙ ЛИТЕРАТУРЫ
• ВВЕДЕНИЕ

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

1. ЧИСЛЕННОЕ РЕШЕНИЕ МЕТОДОМ ЭЙЛЕРА

1.1 Постановка задачи

Решить численно указанную задачу Коши для уравнения первого порядка методом Эйлера, используя табличный процессор Excel.

Исследовать поведение решения на отрезке [0,1.5] с начальным условием y (0)=0, числом отрезков разбиения n=15. Параметр c=0.2.

1.2 Математическая модель задачи (метод Эйлера).

По условию выполнено соотношение:

Пусть начальное значение искомой функции y (x0) =y0. Можно приближенно вычислить следующие значения, находя приращение функции через дифференциал:

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

где i = 0,1,…, n-1, x0, y0 определено из начальных условий, а f (xi, yi) — функция правой части уравнения, вычисленная в узловой точке.

1.3 Исходные данные

Отрезок изменения аргумента: [0,1.5].

Начальное условие y (a): y (0)=0.

Число отрезков разбиения: 15. Тогда шаг h=0.1.

Параметр c=0.2.

1.4 Численное решение уравнения методом Эйлера в Excel

Рассмотрим расчетную таблицу в Excel, содержащую три столбца для значений. Дадим им заголовки x, y,, расположив их в ячейках A2: C2. Для постоянных величин h, c отведем отдельные ячейки E2 и E3. Их заголовки помещены в ячейки D2 и D3.

Расчеты в таблице Excel выполняются по следующему алгоритму:

1. Вычисление первого столбца: первые два значения x = x0 и x1 = x0+h вводятся в ячейки A3 и A4, затем, выделив две эти ячейки, заполняем столбец значений x до достижения конечного значения x=1,5.

2. Затем заполним первую строку расчетной таблицы: в столбце y введем y0 в ячейку B3, в столбце введем в ячейку C3 формулу: =$E$ 2*(A3²-B3²) (вычисляется приращение функции y для текущего значения x в соответствии с формулой).

3. Вычисление второго столбца:

вводим формулу =В3+С3 в ячейку B4 и копируем ее в ячейки B5: B18 (вычисляется новое значение функции y при изменении x на один шаг с помощью линейного приращения по формуле).

4. Вычисление третьего столбца:

копируем формулу из ячейки C3 в ячейки C4: C18

Рис. 1.1. Фрагмент рабочего листа с решением примера по методу Эйлера

2. ЧИСЛЕННОЕ РЕШЕНИЕ МЕТОДОМ РУНГЕ-КУТТА В EXCEL И TURBO PASCAL 7.0

2.1 Постановка задачи

Решить численно указанную задачу Коши для уравнения первого порядка методом Рунге-Кутта, используя табличный процессор Excel и с помощью программы на языке Турбо Паскаль 7.0. Построить графики решений в Excel, с помощью Мастера диаграмм. Провести анализ полученных результатов.

Исследовать поведение решения на отрезке [0,1.5] с начальным условием y (0)=0, числом отрезков разбиения n=15. Параметр c=0.2.

2.2 Математическая модель задачи (метод Рунге-Кутта).

Приводим уравнение к виду .

По условию выполнено соотношение: .

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

где i = 0,1,…, n-1.

2.3 Исходные данные.

Отрезок изменения аргумента: [0,1.5].

Начальное условие y (a): y (0)=1.

Число отрезков разбиения: 15. Тогда шаг h=0.1.

2.4 Численное решение уравнения методом Рунге-Кутта в Excel

Заполняем таблицу Excel в следующем порядке:

1. Первая строка заполнена именами переменных.

2. Ячейка H21 отводится под значения константы h, ячейка G21 — под ее имя.

3. Первый столбец заполняется значениями x.

4. В ячейку B22 вводим значение y0.

5. В ячейки С22, D22, E22, F22 вводим соответственно формулы:

=$H$ 21*(A22²-B22²)

=$H$ 21*((A22+$H$ 21/2)^2-(B22+C22/2)^2)

=$H$ 21*((A22+$H$ 21/2)^2-(B22+D22/2)^2)

=$H$ 21*((A22+$H$ 21)^2-(B22+E22)^2)

6. В ячейку B23 вводим формулу:

=B22+(C22+2*D22+2*E22+F22)/6

7. В ячейки B24: B37 копируем формулу из ячейки B23.

8. В ячейки C22: F36 копируем формулы из ячеек C22: F22.

Рис. 2.1. Фрагмент рабочего листа с решением примера по методу Рунге-Кутта

2.5 Блок-схема алгоритма

Рис. 2.2. Блок-схема алгоритма Рунге-Кутта для дифференциального уравнения 1-го порядка Рис. 2.3. Функция f (x, y: real):real

2.6.Программа на языке Turbo Pascal

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

program rungekutt1_v3;

var

x0,xn, y0, h, x, y, k1, k2,k3,k4:real;

n, i: integer;

inp, ou: text;

function f (x, y: real):real;

begin

f:=exp (x+c*y)+1

end;

begin

assign (inp,'inp1.txt');

reset (inp);

read (inp, x0, xn, n, y0);

close (inp);

h:=(xn-x0)/n;

x:=x0;

y:=y0;

assign (ou,'ou1.txt');

rewrite (ou);

writeln (ou,' x',' ','y');

writeln (ou, x:2:1,' ', y:7:6);

for i:=0 to n-1 do

begin

k1:=h*f (x, y);

k2:=h*f (x+h/2,y+k½);

k3:=h*f (x+h/2,y+k2/2);

k4:=h*f (x+h, y+k3);

y:=y+(k1+2*k2+2*k3+k4)/6;

x:=x+h;

writeln (ou, x:2:1,' ', y:7:6)

end;

close (ou)

end.

2.7 Выполнение расчетов

Исходные значения должны находиться в файле inp1. txt в следующем порядке: левая граница промежутка, правая граница промежутка, число отрезков разбиения, начальное значение функции, параметр.

2.8 Результаты расчетов

Результаты расчетов пишутся в файл ou1.txt.

x y

0.0 0.0

0.1 0.207 392

0.2 0.431 310

0.3 0.674 812

0.4 0.941 709

0.5 1.236 818

0.6 1.566 340

0.7 1.938 427

0.8 2.364 085

0.9 2.858 650

1.0 3.444 396

1.1 4.155 418

1.2 5.47 809

1.3 6.223 828

1.4 7.902 422

1.5 10.714 022

2.9 Представление результатов в виде графиков

Построим графики численного решения данного уравнения в Excel методами Эйлера и Рунге-Кутта.

Рис. 2.4. График зависимости у (х).

2.10 Анализ результатов

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

Результаты, полученные при расчете с использованием Turbo Pascal 7.0 совпадают с расчетами (методами Эйлера и Рунге-Кутта) в Excel, что является косвенным признаком правильности решения задачи.

3. МЕТОД РУНГЕ-КУТТА ДЛЯ ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ ВТОРОГО ПОРЯДКА

3.1 Постановка задачи

1) Решить численно указанную задачу Коши для уравнения второго порядка методом Рунге-Кутта, с помощью программы на языке Turbo Pascal 7.0.

2) Построить графики решения и его производной в Excel с помощью мастера диаграмм.

3) Провести анализ полученных результатов.

3.2 Математическая модель задачи

Рассмотрим уравнение второго порядка, разрешенное относительно второй производной:

(1.12)

на отрезке [a, b] с начальными условиями .

Это уравнение легко свести к системе уравнений первого порядка с помощью замены переменных:. Тогда и уравнение (1.12) сводится к системе первого порядка

(1.13)

с начальными условиями y (a)=y0, z (a)=z0, где за z0 обозначено y10.

Данная система может решаться как система двух уравнений первого порядка (1.8)-(1.9), где

(1.8)

на отрезке [a, b] с начальными условиями

(1.9)

Формулы метода Рунге — Кутта для системы двух уравнений имеют вид:

(1.10)

Где

(1.11)

Последовательно вычисляются на каждом шаге, а затем, где .

3.3 Исходные данные

Уравнение: .

Аргумент: [0.5,1.5].

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

Число отрезков разбиения n: 10.

3.4 Расчет контрольного варианта в Excel для тестирования

Рис. 3.1. Расчет контрольного варианта в режиме отображения данных.

3.5 Блок-схема алгоритма

Рис. 3.2. Блок-схема алгоритма Рунге-Кутта для дифференциального уравнения 2-го порядка.

Рис. 3.3 Функция f1(z:real):real.

Рис. 3.4. Функция f2(x, y, z:real):real.

3.6 Численное решение задачи с использованием Turbo Pascal 7.0

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

program rungekutt2_v3;

var

x0,xn, y0, z0,h, x, y, z, k1,k2,k3,k4,l1,l2,l3,l4,c:real;

n, i: integer;

inp, ou: text;

function f1(z:real):real;

begin

f1:=z

end;

function f2(x, y, z:real):real;

begin

f2:=x-x*z-sin (y)+1

end;

begin

assign (inp,'inp2.txt');

reset (inp);

read (inp, x0, xn, n, y0,z0);

close (inp);

h:=(xn-x0)/n;

x:=x0;

y:=y0;

z:=z0;

assign (ou,'ou2.txt');

rewrite (ou);

writeln (ou,' x',' ','y',' ','z');

writeln (ou, x:2:1,' ', y:7:6,' ', z:7:6);

for i:=0 to n-1 do

begin

k1:=h*f1(z);

l1:=h*f2(x, y, z);

k2:=h*f1(z+l½);

l2:=h*f2(x+h/2,y+k½, z+l½);

k3:=h*f1(z+l2/2);

l3:=h*f2(x+h/2,y+k2/2,z+l2/2);

k4:=h*f1(z+l3);

l4:=h*f2(x+h, y+k3,z+l3);

y:=y+(k1+2*k2+2*k3+k4)/6;

z:=z+(l1+2*l2+2*l3+l4)/6;

x:=x+h;

writeln (ou, x:2:1,' ', y:7:6,' ', z:7:6)

end;

close (ou)

end.

3.7 Выполнение расчетов

Исходные значения должны находиться в файле inp2. txt в следующем порядке:

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

3.8 Результаты расчетов

Результаты расчетов пишутся в файл ou2.txt.

x y z

0.5 1.0 1.0

0.6 1.100 693 1.12 891

0.7 1.202 401 1.20 476

0.8 1.304 646 1.23 823

0.9 1.407 064 1.24 112

1.0 1.509 408 1.22 562

1.1 1.611 554 1.20 359

1.2 1.713 494 1.18 593

1.3 1.815 320 1.18 223

1.4 1.917 212 1.20 041

1.5 2.19 422 1.24 669

3.9 Представление результатов в виде графиков

Построим графики численного решения данного уравнения в Excel методом Рунге-Кутта.

Рис. 3.5. График решения уравнений системы.

3.10 Анализ результатов

Результаты, полученные при расчете с использованием Turbo Pascal 7.0, совпадают с контрольным вариантом расчета по методу Рунге-Кутта, произведенном в Excel, что является косвенным признаком правильности решения.

4. МОДЕЛЬ ТИПА «ХИЩНИК-ЖЕРТВА» С УЧЕТОМ ВНУТРИВИДОВОГО ВЗАИМОДЕЙСТВИЯ

4.1 Постановка задачи

1) Решить численно задачу Коши для системы дифференциальных уравнений первого порядка, являющейся моделью «хищник — жертва» при различных начальных условиях и различных значениях параметров системы. Величину интервала решения и шаг подобрать самостоятельно, так чтобы решение не выходило за пределы сотен единиц и давало картину поведения траекторий. Решить задачу на отрезке [0,T] .

2) Решение выполнить с помощью системы Турбо Паскаль 7.* или Delphi.

3) Контрольный вариант рассчитать самостоятельно вручную или с помощью Excel.

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

5) Провести анализ поведения решений.

4.2 Математическая модель задачи

Система уравнений, описывающая динамику популяций двух видов, взаимодействующих между собой по типу «хищник-жертва» и с учетом внутривидового взаимодействия:

Заданные начальные условия при t=0

4.3 Исходные данные

Начальное и конечное значение аргумента t: t0, tn.

Начальные значения неизвестных функций .

Число отрезков разбиения области интегрирования n.

t0=0, tn= 15, n=30, x0=100, x0=90, x0=85, y0=30.

Параметры:

a=0.25, b=-0.01, c=0.01, d=-1.0

Начальные данные: (100,30), (90,30), (85,30).

Таким образом, получаем систему:

4.4 Расчет контрольного варианта в Excel для тестирования.

Рис. 4.1. Расчет контрольного варианта для тестирования (для демонстрации приведено только 3 точки; специальная вставка).

4.5 Блок-схема алгоритма

Рис. 4.2. Блок-схема алгоритма решения уравнения второго порядка.

Рис. 4.3. Функция f1(x, y: real):real

Рис. 4.4. Функция f2(x, y: real):real

4.6 Численное решение задачи с использованием Turbo Pascal 7.0

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

program rungekutt3_v3;

var

t0,tn, x0, xi, y0, yi, h, t, x, y, k1, k2,k3,k4,l1,l2,l3,l4,a, b, c, d: real;

n, i: integer;

inp, ou: text;

function f1(x, y: real):real;

begin

f1:=a*x+b*x*y

end;

function f2(x, y: real):real;

begin

f2:=d*y+c*x*y

end;

begin

assign (inp,'inp3.txt');

reset (inp);

read (inp, t0, tn, n, x0,y0,a, b, c, d);

close (inp);

h:=(tn-t0)/n;

t:=t0;

x:=x0;

y:=y0;

assign (ou,'ou3.txt');

rewrite (ou);

writeln (ou,' t',' ','x',' y');

writeln (ou, t:2:1,' ', x:7:6,' ', y:7:6);

for i:=0 to n-1 do

begin

k1:=h*f1(x, y);

l1:=h*f2(x, y);

k2:=h*f1(x+k½, y+l½);

l2:=h*f2(x+k½, y+l½);

k3:=h*f1(x+k2/2,y+l2/2);

l3:=h*f2(x+k2/2,y+l2/2);

k4:=h*f1(x+k3,y+l3);

l4:=h*f2(x+k3,y+l3);

x:=x+(k1+2*k2+2*k3+k4)/6;

y:=y+(l1+2*l2+2*l3+l4)/6;

t:=t+h;

writeln (ou, t:2:1,' ', x:7:6,' ', y:7:6)

end;

close (ou)

end.

4.7 Выполнение расчетов

Исходные значения должны находиться в файле inp3. txt в следующем порядке:

t0, tn, n, x0, y0, a, b, c, d.

4.8 Результаты расчетов

Результаты расчетов пишутся в файл ou3.txt.

t x y

0.0 100.0 30.0

0.5 97.561 239 29.815 776

1.0 95.353 243 29.288 566

1.5 93.509 932 28.479 413

2.0 92.123 606 27.469 108

2.5 91.246 723 26.344 787

3.0 90.897 734 25.189 133

3.5 91.68 114 24.73 624

4.0 91.728 650 23.55 723

4.5 92.834 069 22.179 046

5.0 94.325 805 21.475 262

5.5 96.133 177 20.966 690

6.0 98.173 439 20.668 779

6.5 100.351 315 20.591 967

7.0 102.558 763 20.742 577

7.5 104.675 808 21.122 501

8.0 106.573 439 21.727 539

8.5 108.119 585 22.544 379

9.0 109.188 919 23.546 528

9.5 109.676 451 24.690 050

10.0 109.513 570 25.910 610

10.5 108.683 393 27.123 941

11.0 107.230 929 28.231 615

11.5 105.263 667 29.132 734

12.0 102.940 443 29.739 682

12.5 100.450 388 29.993 707

13.0 97.987 412 29.875 306

13.5 95.726 963 29.406 092

14.0 93.810 121 28.642 112

14.5 92.336 705 27.661 622

15.0 91.365 984 26.551 514

t x y

0.0 90.0 30.0

0.5 88.124 156 28.394 467

1.0 87.12 742 26.674 766

1.5 86.658 548 24.967 828

2.0 87.24 499 23.371 275

2.5 88.55 005 21.953 894

3.0 89.682 440 20.760 590

3.5 91.829 580 19.818 959

4.0 94.409 030 19.145 673

4.5 97.320 635 18.751 745

5.0 100.447 745 18.646 266

5.5 103.653 281 18.838 414

6.0 106.776 833 19.337 464

6.5 109.634 518 20.150 400

7.0 112.23 954 21.276 600

7.5 113.737 034 22.699 264

8.0 114.582 555 24.374 201

8.5 114.418 204 26.218 489

9.0 113.186 582 28.104 125

9.5 110.944 203 29.863 405

10.0 107.869 782 31.310 759

10.5 104.242 954 32.278 633

11.0 100.396 886 32.655 897

11.5 96.660 823 32.413 057

12.0 93.312 336 31.604 059

12.5 90.551 709 30.346 190

13.0 88.499 048 28.788 864

13.5 87.206 469 27.83 737

14.0 86.676 001 25.363 928

14.5 86.876 414 23.734 166

15.0 87.755 735 22.269 710

t x y

0.0 85.0 30.0

0.5 83.375 891 27.708 055

1.0 82.724 512 25.446 600

1.5 82.981 548 23.347 415

2.0 84.68 146 21.494 220

2.5 85.900 092 19.933 744

3.0 88.390 379 18.688 092

3.5 91.447 368 17.765 652

4.0 94.970 044 17.169 616

4.5 98.841 344 16.904 054

5.0 102.920 315 16.977 785

5.5 107.34 024 17.406 067

6.0 110.970 899 18.209 699

6.5 114.478 409 19.410 485

7.0 117.269 715 21.21 441

7.5 119.45 257 23.30 314

8.0 119.534 328 25.376 914

8.5 118.555 694 27.929 736

9.0 116.83 872 30.474 461

9.5 112.293 600 32.730 777

10.0 107.552 495 34.405 945

10.5 102.352 137 35.271 435

11.0 97.202 394 35.227 990

11.5 92.534 733 34.326 265

12.0 88.649 307 32.736 884

12.5 85.710 787 30.691 654

13.0 83.775 361 28.425 558

13.5 82.827 091 26.138 116

14.0 82.809 436 23.977 472

14.5 83.646 711 22.41 294

15.0 85.255 835 20.386 485

4.9 Представление результатов в виде графиков.

Построим графики в Excel численного решения данного уравнения.

Рис. 4.1. Траектория решения при различных исходных данных.

4.10 Анализ результатов

Из графического представления решения системы дифференциальных уравнений первого порядка, описывающей динамику популяций двух видов, взаимодействующих между собой по типу «хищник-жертва» и с учетом внутривидового взаимодействия, видно, что поведение системы изменяется при разных параметрах. В данном случае при изменении y (начального количества хищников) изменяется амплитуда колебаний численности обоих популяций в системе. Тем не менее, во всех трех случаях система ведет себя стабильно.

Результаты, полученные при расчете с использованием Turbo Pascal 7.0, совпадают с контрольным вариантом, рассчитанном в Excel, что является косвенным признаком правильности решения.

ЗАКЛЮЧЕНИЕ

В данной работе рассмотрены методы решения дифференциальных уравнений и систем. Закреплены навыки решения дифференциальных уравнений и систем первого порядка методом Рунге-Кутта, которые могут быть применены в будущем при построении различных математических моделей и в прикладных задачах. В ходе решения использовались программы: Turbo Pascal 7.1. процессор MS Excel. Отчет оформлен в Microsoft Word.

численный метод дифференциальное уравнение

Н.С. Пискунов. Дифференциальное и интегральное исчисление.

Лабораторные работы по курсу «Вычислительная математика и применение ЭВМ», методическое пособие. — Ленинград, 1987. — 160 с.

Турчак Л. И. Основы численных методов. — М.: «Наука», 1987. ;

Воробьева Г. Н., Данилова А. Н. Практикум по вычислительной математике, М.:"Высшая школа", 1991. 208 с.

Информатика. Программирование в среде Турбо Паскаль 7.0. Лабораторные работы 1−3, 4−6, 7−9. СПб.: СПГГИ, 2003.

Турбо Паскаль 7.0. Киев, Торгово-издательское бюро BHV, 1997.

Вычислительная техника и программирование. Практикум по программированию. — М.:"Высшая школа", 1990 / Под ред. А. В. Петрова — 400 с.