Лабораторная работа №4 численное интегрирование. ЧИсленное решение обыкно-венных дифференциальных уравнений

Численное интегрирование

Если функция f(x) непрерывна на [a,b] и известна ее первообразная Ф(х), то определенный интеграл от этой функции вычисляется по формуле Ньютона-Лейбница:

. (10)

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

Обычный прием численного интегрирования состоит в том, что данную функцию f(x) на рассматриваемом промежутке заменяют интерполирующей функцией простого вида F(x), а затем приближенно полагают:

. (11)

Метод прямоугольников

Разобьем отрезок [a,b] на n равных промежутков точками x₀, x₁,…, x_n. Величина h = (b-a)/n – длина каждого промежутка разбиения – называется шагом интегрирования.

Заменим функцию f(x) на каждом промежутке постоянной функцией, принимающей значение, равное значению в левом (правом) конце промежутка. Получим при этом формулу левых (правых) прямоугольников (как площадь ступенчатой фигуры, образованной прямоугольниками).

Формула левых прямоугольников:

. (12)

Формула правых прямоугольников:

. (13)

Справедлива следующая оценка погрешности формул прямоугольников:

;

(14)

Формула трапеции

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

Формула трапеции имеет вид:

. (15)

Оценка точности формулы:

;

. (16)

Формула Симпсона

Формула Симпсона для многочленов 2 степени имеет вид:

(17)

Формула остаточного члена:

;

(18)

Задание 4.1. Вычислить значение определенного интеграла, заданного в таблице 2.

Таблица 2 – Задания для расчетов

Вариант	Формула интеграла		Вариант	Формула интеграла
1	1)		7	1)
	2)			2)
2	1)		8	1)
	2)			2)

Продолжение таблицы

3	1)		9	1)
	2)			2)
4	1)		10	1)
	2)			2)
5	1)		11	1)
	2)			2)
6	1)		12	1)
	2)			2)

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

Дифференциальные уравнения и их системы можно решать численными методами (в частности, методом конечных разностей).

Решение заключается в нахождении ряда значений x_i и y_i искомой зависимости y(x) при i, изменяющемся от 0 до N при шаге изменения х, равном h. Будем рассматривать способы решения дифференциального уравнения, при которых h = const.

Рассмотрим методы решения обыкновенных дифференциальных уравнений первого порядка вида Y' = f(X,Y).

Mетод Эйлера

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

x_i+1=x₁+h; (27)

y_i+1=y₁+h*f(x,y).

При этом для i = 0 значения x₀ и y₀ должны быть известны как начальные условия. Погрешность метода пропорциональна h².

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

x_i+1=x₁+h;

y_i+1=y₁+h*(f(x₁,y₁)+f(x₁+h,y₁+h*f(x₁,y₁)))/2 (28)

Погрешность метода пропорциональна h³.

Метод Рунге-Кутта

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

k1 = h* f(x_i,y_i);

k2 = h* f(x_i+h/2,y_i+k1/2);

k3 = h* f(x_i+h/2,y_i+k2/2); (29)

k4=h*f(x_i+h,y_i+k3) ;

x_i+1 = x_i + h;

y_i+1 = y_i + (k1+2*k2+2*k3+k4)/6.

Погрешность метода пропорциональна h⁵.

Для решения обыкновенных дифференциальных уравнений 2-го порядка

Y''+p(x)Y'+g(x)Y=f(x) с граничными условиями:

k11*Y(a)+k12*Y' (a)=A,

k21*Y(b)+k22*Y´(b)=B

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

; ;; (30)

. (31)

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

ПРИМЕР. Решить обыкновенное дифференциальное уравнение 2-го порядка:

y''+xy'-0.5y/x=1

Решение:

Разбив отрезок [2, 2.3] на части с шагом h = 0.1, получим четыре узловые точки с абсциссами x₀=2, x₁=2.1, x₂=2.2, x₃=2.3. Две точки x₀=2 и x₃=2.3 являются конечными, две другие – внутренними. Данное уравнение во внутренних точках заменим конечно-разностным:

(i=1, 2).

Из краевых условий составим конечно-разностные уравнения в конечных точках:

,

y3=2.15.

Задача сводится к решению системы уравнений:

Примеры решения дифференциальных уравнений с помощью встроенных функций системы MathСad приведены в приложениях Д, Е.

Задание 7.1. Решите на отрезке [x₀,x_end] задачу Коши методом Рунге-Кутта с постоянным шагом. Вид уравнения и начальные значения заданы в таблице 3. Изобразите графики решений, вычисленных с шагамиh , 2h и h/2.

Таблица 3 – Данные для расчета

Вариант	Функция	x0	xend
1		2	3	1	0.1
2		3	4	1	0.1
3		0	1	2	0.1
4		2	3	1	0.1
5		1	2	1	0.1
6		0	1	1	0.1
7		0	1	2	0.1
8		0	1	1	0.1
9		2	3	2	0.1
10		0	1	3	0.1
11		0	2	1	0.1
12		1	3	1	0.1

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

1. Присвойте начальное значение решения переменной у₀ .

2. Определите правую часть уравнения f(x,y).

3. Вычислите решение, используя функцию rkfixed с параметром N, вычисленным по формуле N= .

4. Сохраните решение в матрице У1.

5. Вычислите решение, используя функцию rkfixed с параметром N, вычисленным по формуле N= 2.

6. Сохраните решение в матрице У2.

7. Вычислите решение, используя функцию rkfixed с параметром N, вычисленным по формуле N= .

8. Сохраните решение в матрице У3.

9. Постройте на одном графике все три найденных решения.

10. Оцените погрешности найденных решений по формуле Рунге.

Задание 7.2. Решите задачу Коши

на отрезке [a,b] методом Рунге-Кутты с постоянным шагом h=0.1. Изобразите графики решений, вычисленных с шагом h, 2h и h/2. Вид уравнений и начальные значения заданы в таблице 4.

Таблица 4 – Данные для расчета

Вариант					a	b
1			1	0	-1	1
2			0.5	1.5	0	2
3			-1	1	0	4
4			1	0	0	5

Продолжение таблицы

5			0.2	0	-1	1
6			0	0	0	4
7			0.5	-0.5	-1	3
8			-0.6	2	2	5
9			0	0	-1	3
10			1.2	1.2	0	2
11			1	1	1	3
12			0.8	3.5	2	4

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

1. Присвойте переменной ORIGIN значение, равное единице.

2. Присвойте начальное значение решения переменной у₀ .

3. Определите правую часть уравнения f(x,y).

4. Вычислите решение, используя функцию rkfixed с параметром N, вычисленным по формуле N= .

5. Сохраните решение в матрице У1.

6. Вычислите решение, используя функцию rkfixed с параметром N, вычисленным по формуле N= 2.

7. Сохраните решение в матрице У2.

8. Вычислите решение, используя функцию rkfixed с параметром N, вычисленным по формуле N= .

9. Сохраните решение в матрице У3.

10. Постройте на одном графике все три найденные решения.

11. Оцените погрешности найденных решений по формуле Рунге.

Задание 7.3. Найдите общее решение линейного однородного уравнения второго порядка . Решите задачу Коши,. Изобразите его график. Значения параметрова₁, а₂ и а заданы в таблице 5.

Таблица 5 – Данные для расчета

Вариант	a1	a2	y(a)	y(a)	a
1	2	0	1	1	0
2	-4	4	0	1	0
3	2	0	0	0	/2
4	0	1	3	0	-/2
5	2	5	0	0	1
6	-4	4	1	1	0.3
7	6	13	-1	1	0.25
8	0	1	4	1	/2
9	2	5	6	2	-/2
10	-4	8	0	2	1
11	2	0	0	2	0
12	-4	4	-1	0.5	2

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

1. Присвойте переменной ORIGIN значение, равное единице.

2. Присвойте начальное значение решения переменной у₀ .

3. Определите правую часть уравнения f(x,y).

4. Вычислите решение, используя функцию rkfixed с параметром N, вычисленным по формуле N= .

5. Сохраните решение в матрице У1.

Задание 7.4. Постройте графики решения и фазовые портреты динамической системы

моделирующей взаимодействие популяций при заданных значениях параметров a,b,c,d . Значения параметров заданы в таблице 6. Исследуйте поведение решения, изменяя параметры.

Таблица 6 – Данные для расчета

Вариант	a	b	c	d
1	4	3	2	1
2	4	3	2	2
3	4	3.5	2	1
4	4	3.5	2	2
5	4	3.5	3	1
6	4	3.5	3	2
7	4	3.5	4	4
8	5	4.5	2	1
9	5	4	2	1
10	5	4	3	1
11	5	4	2	2
12	5	4	2	3

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

1. Присвойте переменной ORIGIN значение, равное единице.

2. Определите вектор-столбец начальных условий для первой задачи Коши.

3. Определите вектор-столбец правых частей системы.

4. Выберите значение шага интегрирования h и вычислите количество шагов N интегрирования системы на отрезке [x₀,x_end] по формуле N= .

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

6. Изобразите соответствующую фазовую кривую и график решения.

7. Определите векторы-столбцы начальных условий для каждого начального условия.

8. Решите соответствующие задачи, сохранив каждое решение в отдельной матрице.

9. Изобразите соответствующие фазовые кривые и графики решений.

Задание 7.5. Исследуйте поведение системы

моделирующей взаимодействие популяций. Выполните вычисления для значений a,b,c,d из задания 7.4 для приведённых ниже значений  (таблица 7).

Таблица 7 – Данные для расчета

Вариант		Вариант		Вариант	
1	0.1	5	0.05	9	0.17
2	0.15	6	0.22	10	0.18
3	0.20	7	0.12	11	0.1
4	0.25	8	0.14	12	0.2