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

Размещено на http://www.allbest.ru/

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

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

СОДЕРЖАНИЕ

• ВВЕДЕНИЕ
• 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^2-B3^2) (вычисляется приращение функции 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^2-B22^2)

=$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+k1/2);

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.000000

0.1 0.207392

0.2 0.431310

0.3 0.674812

0.4 0.941709

0.5 1.236818

0.6 1.566340

0.7 1.938427

0.8 2.364085

0.9 2.858650

1.0 3.444396

1.1 4.155418

1.2 5.047809

1.3 6.223828

1.4 7.902422

1.5 10.714022

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+l1/2);

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

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.000000 1.000000

0.6 1.100693 1.012891

0.7 1.202401 1.020476

0.8 1.304646 1.023823

0.9 1.407064 1.024112

1.0 1.509408 1.022562

1.1 1.611554 1.020359

1.2 1.713494 1.018593

1.3 1.815320 1.018223

1.4 1.917212 1.020041

1.5 2.019422 1.024669

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+k1/2,y+l1/2);

l2:=h*f2(x+k1/2,y+l1/2);

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.000000 30.000000

0.5 97.561239 29.815776

1.0 95.353243 29.288566

1.5 93.509932 28.479413

2.0 92.123606 27.469108

2.5 91.246723 26.344787

3.0 90.897734 25.189133

3.5 91.068114 24.073624

4.0 91.728650 23.055723

4.5 92.834069 22.179046

5.0 94.325805 21.475262

5.5 96.133177 20.966690

6.0 98.173439 20.668779

6.5 100.351315 20.591967

7.0 102.558763 20.742577

7.5 104.675808 21.122501

8.0 106.573439 21.727539

8.5 108.119585 22.544379

9.0 109.188919 23.546528

9.5 109.676451 24.690050

10.0 109.513570 25.910610

10.5 108.683393 27.123941

11.0 107.230929 28.231615

11.5 105.263667 29.132734

12.0 102.940443 29.739682

12.5 100.450388 29.993707

13.0 97.987412 29.875306

13.5 95.726963 29.406092

14.0 93.810121 28.642112

14.5 92.336705 27.661622

15.0 91.365984 26.551514

t x y

0.0 90.000000 30.000000

0.5 88.124156 28.394467

1.0 87.012742 26.674766

1.5 86.658548 24.967828

2.0 87.024499 23.371275

2.5 88.055005 21.953894

3.0 89.682440 20.760590

3.5 91.829580 19.818959

4.0 94.409030 19.145673

4.5 97.320635 18.751745

5.0 100.447745 18.646266

5.5 103.653281 18.838414

6.0 106.776833 19.337464

6.5 109.634518 20.150400

7.0 112.023954 21.276600

7.5 113.737034 22.699264

8.0 114.582555 24.374201

8.5 114.418204 26.218489

9.0 113.186582 28.104125

9.5 110.944203 29.863405

10.0 107.869782 31.310759

10.5 104.242954 32.278633

11.0 100.396886 32.655897

11.5 96.660823 32.413057

12.0 93.312336 31.604059

12.5 90.551709 30.346190

13.0 88.499048 28.788864

13.5 87.206469 27.083737

14.0 86.676001 25.363928

14.5 86.876414 23.734166

15.0 87.755735 22.269710

t x y

0.0 85.000000 30.000000

0.5 83.375891 27.708055

1.0 82.724512 25.446600

1.5 82.981548 23.347415

2.0 84.068146 21.494220

2.5 85.900092 19.933744

3.0 88.390379 18.688092

3.5 91.447368 17.765652

4.0 94.970044 17.169616

4.5 98.841344 16.904054

5.0 102.920315 16.977785

5.5 107.034024 17.406067

6.0 110.970899 18.209699

6.5 114.478409 19.410485

7.0 117.269715 21.021441

7.5 119.045257 23.030314

8.0 119.534328 25.376914

8.5 118.555694 27.929736

9.0 116.083872 30.474461

9.5 112.293600 32.730777

10.0 107.552495 34.405945

10.5 102.352137 35.271435

11.0 97.202394 35.227990

11.5 92.534733 34.326265

12.0 88.649307 32.736884

12.5 85.710787 30.691654

13.0 83.775361 28.425558

13.5 82.827091 26.138116

14.0 82.809436 23.977472

14.5 83.646711 22.041294

15.0 85.255835 20.386485

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 с.

Размещено на Allbest.ru