Mathcad - ЛР6

Министерство образования Российской Федерации Санкт - Петербургский государственный политехнический университет

Институт менеджмента и информационных технологий

Кафедра ПО ВТ и АС

Лабораторный практикум по курсу вычислительной математики

ОТЧЕТ

по лабораторной работе №6

Тема:

ЧИСЛЕННОЕ РЕШЕНИЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ (ЗАДАЧА КОШИ)

Выполнил: _______________

Группа: _____

Вариант: 14

Проверил: Царев В.А.

Отметка о зачете __________________

" ___ " ______________ 2006 г.

Череповец

2006

Задача 6.1. Найти приближенное решение задачи Коши для обыкновенного дифференциального

уравнения (ОДУ) 1 порядка

y '(t) = f (t, y(t)), t [t0,T ], (1) y(t0) = y0

и оценить погрешность решения задачи. ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ:

1.Задать исходные данные: функцию f правой части, начальное значение y0 .

2.Используя функцию eyler (см. ПРИЛОЖЕНИЕ 6.B), найти приближенное решение задачи Коши с шагом h=0.1 по явному методу Эйлера.

3.Используя встроенную функцию rkfixed пакета MATHCAD, найти приближенное решение задачи Коши с шагом h=0.1 (см. ПРИЛОЖЕНИЕ 6.B).

4.Найти решение задачи Коши аналитически.

5.Построить таблицы значений приближенных и точного решений. На одном чертеже построить графики приближенных и точного решений.

6.Оценить погрешность приближенных решений двумя способами:

a) по формуле ε = max | y(ti ) − yi |; здесь y(ti )и yi - значения точного и приближенного

0≤i≤N

решений в узлах сетки ti , i=1,..N;

b) по правилу Рунге (по правилу двойного пересчета) (см. ПРИЛОЖЕНИЕ 6.C).

7.Выяснить, при каком значении шага h=h* решение, полученное по методу Эйлера, будет иметь такую же погрешность (см. п. 6а), как решение, полученное с помощью метода Рунге-Кутты с шагом h=0.1.

УКАЗАНИЕ. В п. 7 рекомендуется провести серию вычислений решения по методу Эйлера, дробя шаг h пополам.

Теоретический материал:

Метод Эйлера относится к явным методам и реализуется следующей формулой: yk +1 = yk + h × f (xk , yk )

где f (xk , yk ) - производная функции.

Метод Эйлера имеет второй порядок локальной точности и первый порядок глобальной точности.

Значительно более точный результат численного дифференцирования дает метод Рунге- Кутты четвертого порядка точности. Вычисления по данному методу производятся по следующей схеме:

yi+1	= yi +		h	(K1 + 2K2 + 2K3 + K4 )
		6
K1 = f (xi , yi )
K2	= f (xi	+ 0.5h, yi			+ 0.5h × K1 )
K3	= f (xi	+ 0.5h, yi			+ 0.5h ×K 2)
K4	= f (xi + h, yi + h × K3 )

Правило Рунге практической оценки погрешности (правило двойного пересчета):

Одним из наиболее простых, широко применяемых и достаточно эффективных методов

оценки погрешности и уточнения полученного решения при численном решении ОДУ является правило Рунге. Для оценки погрешности решения дифференциального уравнения (разницы между точным и приближенным решением в узлах) применяют первую формулу Рунге:

	y(h / 2)	− y(h)						y(ti ) − yi(h / 2) ≈ εi(h) ,
где εi(h) =
	2i	i		, i= 0, … , n – задает номера узлов сетки ti , построенной с шагом h ,
	2 p −1
			yi(h)		и y2(hi / 2) - приближенные решения уравнения с шагом h и h / 2
p – порядок метода,
соответственно в узле ti .
Исходные данные:
Правая часть уравнения:					f(t , y) :=	y		-	12
						t
									t	3

Начальное значение:

Решение задачи:

y0 := 4

Отрезок: t0 := 1

T := 2 Шаг сетки: h := 0.1

Число узлов сетки: N :=	T - t0	N = 10

h

i := 1 .. N

ti := t0 + i×h

Функция, реализующая явный метод Эйлера; возвращает вектор решения:

eyler(f , y0, t0 , h, N) := y0 ¬ y0

for i Î 0 .. N - 1

yi+1 ¬ yi + h×f(t0 + i×h, yi)

y

yE := eyler(f , y0 , t0 , h, N)

Входные параметры:

f - функция правой части; y0 - начальное значение; t0 - начальная точка отрезка;

h - шаг сетки;

N - число узлов сетки.

			0
		0	4
		1	3.2
		2	2.5893313298272
		3	2.11066449620169
	yE =	4	1.72682394083441
		5	1.4128507237803
		6	1.15148521647677
		7	0.930484292506567
		8	0.740968713539418
		9	0.576372325196958
		10	0.431755093734077

ВЫЧИСЛЕНИЕ РЕШЕНИЯ ПРИ ПОМОЩИ ВСТРОЕННОЙ ФУНКЦИИ (она реализует метод Рунге-Кутты 4-го порядка):

yRK4 := rkfixed(y, t0 , T , N , f)	- встроенная функция;
	входные параметры:
	y - вектор начальных значений;
	t0- начальная точка отрезка;
	T - конечная точка отрезка;
	N - число узлов сетки;
	f - функция правой части.

Функция rkfixed возвращает матрицу, первый столбец которой содержит узлы сетки, а второй - приближенное решение в этих узлах.

		0		1
	0		1	4
	1		1.1	3.30576298486409
	2		1.2	2.77774195868259
	3		1.3	2.36681861272695
yRK4 =	4		1.4	2.04076377359122
	5		1.5	1.77771918609449
	6		1.6	1.56243607225608
	7		1.7	1.38401420013955
	8		1.8	1.23449439708241
	9		1.9	1.10795524063915
	10		2	0.999917609225951

	yTrue(t) :=	4		- точное решение диф. уравнения
		t2

	3
yEi
yRK4i , 1	2
yTrue(ti)
	1
	1	1.2	1.4	1.6	1.8
				ti
Произведем оценку погрешностей приближенных решений:

Вычисление погрешностей решения:

eE_1i :=		yTrue(ti) - yEi			eRK4_1i :=		yTrue(ti) - yRK4i , 1
max_e_E := max(eE_1)					max_e_RK4 := max(eRK4_1)
max_e_E = 0.568244906265923					max_e_RK4 = 8.23907740494789 ´ 10- 5

Вычисление погрешностей решения по правилу Рунге:

Вычисление приближенных решений с шагом h/2:

h2 :=

h

N2 :=

T - t0

N2 = 20

2

h2

yEh2 :=

eyler(f , y0 , t0 , h2 , N2)

yRK4h2 := rkfixed(y, t0 , T , N2 , f)

Вычисление погрешностей по точкам:

i := 0 .. N

zEi :=

yEh22×i - yEi

zRK4i :=

(yRK4h2á1ñ)

2×i

- (yRK4á1ñ )

i

15

Значение погрешностей:

max(zE) = 0.285306411369347

max(zRK4) = 5.14246690531313 ´ 10- 6

Расчет абсолютной погрешности через функцию Mmax:

Mmax(y, yr , n) := max ¬ 0 for i Î 0 .. n

max ¬ yi - yri if yi - yri > max

max

Функция, рассчитывающая текущую погрешность и соответствующий ей шаг при заданной целевой погрешности.

get_h(h , e , yTrue) :=	s0 ¬ 1
	while s0 > e
		N ¬			T - t0
						h
		for		i Î 0 .. N
			ti ¬ t0 + i×h
			yti ¬ yTrue(ti)
		yE	¬ eyler(f , y0 , t0 , h, N)
		s0 ¬ Mmax(yt , yE , N)
		s1 ¬ h
		h ¬		h
				2
	s

æ	6.86646663453638 ´	− 5	öпогрешность
ç		10	÷
get_h(0.1 , max_e_RK4 , yTrue) = ç	1.220703125 ´ 10− 5		÷
è			ø
			шаг

При шаге 1.220703125 ´ 10− 5 погрешность метода Эйлера станет лучше погрешности

метода Рунге-Кутты четвертого порядка с шагом h=0.1

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

Задача 6.2. Задача Коши для ОДУ 2 порядка

mx′′ + Hx′ + kx = f (t) , x(0) = x0

x′(0) = v0

t [0,T ],

описывает движение груза массы m, подвешенного к концу пружины. Здесь x(t) – смещение груза от положения равновесия, H – константа, характеризующая силу сопротивления среды, k –коэффициент

упругости пружины, f(t) – внешняя сила. Начальные условия: x0 – смещение груза в начальный момент

времени t=0, v0 – скорость груза в начальный момент времени. Промоделировать движение груза на

временном отрезке [0,T] при заданных в индивидуальном варианте трех наборах (I, II, III) значений параметров задачи. Для каждого набора по найденной таблице (или графику) решения задачи определить максимальное и минимальное значения функции x(t) и моменты времени, в которые эти значения достигаются. Предложить свой вариант задания параметров, при которых характер колебаний груза существенно отличается от рассмотренного ранее.

ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ:

1. Заменить исходную задачу эквивалентной задачей Коши для системы ОДУ 1 порядка:

	′	= x2
	x1
	′	=	f (t) − Hx2 − kx1
	x2
			m	(2)

x1(0) = x0 x2 (0) = v0

2.Для каждого варианта выбора параметров решить задачу (2) с помощью встроенной функции rkfixed пакета MATHCAD (см. ПРИЛОЖЕНИЕ 7.B) с шагом h=0.1.

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

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

Исходные данные I:

	H (t) := 1	K := 1 m := 1		f(t) := cos(t)
	X0 := 0	V0 := 0	t0 := 0	T := 34 Шаг сетки: h := 0.1
	Решение задачи I:
	Число узлов сетки:		N :=	T - t0
				h

Формирование вектора правой части системы ОДУ (1) и вектора начальных условий для применения встроенной функции rkfixed:

f1(t , x1 , x2) := x2			f2(t , x1 , x2) :=			f(t) - H(t) ×x2 - K×x1
						m
	æX0 ö		æf1(t , x0 , x1) ö
x :=	çV0	÷	D (t , x) := çf2	(t , x , x ) ÷
	è	ø	è	0 1	ø

xRK4 := rkfixed(x, t0 , T , N , D)

get_minimum (m , n1 , n2) :=			min ← 1015
			for	i n1 .. n2
			if	mi , 1 < min
				min ← mi , 1
				result0 ← mi , 0
				result1 ← mi , 1
				result2 ← i
			result
get_maximum(m , n1 , n2) :=			max ← 0
			for	i n1 .. n2
			if	mi, 1 > max
				max ← mi , 1
				result0 ← mi, 0
				result1 ← mi, 1
				result2 ← i
			result
График решения:			i := 0 .. N
	1
xRK4i , 1	0
	1 0	10	20	30
		xRK4i , 0
Исходные данные II:

æ	17.3	ö
get_minimum (xRK4, 0 , 340) = ç	−0.999910402134673	÷
ç		÷
è	173	ø

æ26.7000000000001 ö get_maximum(xRK4, 0 , 340) = çç0.99999715251772 ÷÷ è 267 ø

	H (t) := 1	K := 1 m := 14		f(t) := cos(t)
	X0 := 0	V0 := 0	t0 := 0	T := 34 Шаг сетки: h := 0.1
	Решение задачи II:
	Число узлов сетки:		N :=	T − t0
				h

Формирование вектора правой части системы ОДУ (1) и вектора начальных условий для применения встроенной функции rkfixed:

	f1(t , x1 , x2) := x2	f2(t , x1 , x2) :=	f(t) - H(t) ×x2 - K×x1
			m

æX0 ö			æf1	(t , x , x ) ö
		D (t , x) := çf2		0	1
x := çV0	÷			(t , x , x ) ÷
è	ø		è	0	1 ø
xRK4 :=	rkfixed(x, t0 , T , N , D)
График решения:			i := 0 .. N
		0.1
xRK4i , 1		0
		0.1
		0	20
			xRK4i , 0
Исходные данные III:

æ	12.4	ö
get_minimum (xRK4, 0 , 340) = ç-0.123912023170379 ÷
ç		÷
è	124	ø
æ	2.9	ö
get_maximum(xRK4, 0 , 340) = ç	0.117792615900578	÷
ç		÷
è	29	ø

	H (t) := 1	K := 1 m := 28		f(t) := cos(t)
	X0 := 0	V0 := 0	t0 := 0	T := 34 Шаг сетки: h := 0.1
	Решение задачи III:
	Число узлов сетки:		N :=	T - t0
				h

Формирование вектора правой части системы ОДУ (1) и вектора начальных условий для применения встроенной функции rkfixed:

f1(t , x1 , x2) := x2			f2(t , x1 , x2) :=				f(t) - H(t) ×x2 - K×x1
							m
	æX0	ö	æf1	(t , x , x ) ö
			D (t , x) := çf2	0	1
x :=	çV0	÷		(t , x , x ) ÷
	è	ø	è	0	1	ø

xRK4 := rkfixed(x, t0 , T , N , D)

График решения:		i := 0 .. N
xRK4i , 1	0
	0	20
		xRK4i , 0
Исходные данные IV:

æ	18.7	ö
get_minimum (xRK4, 0 , 340) = ç-0.060414912262805 ÷
ç		÷
è	187	ø
æ	3	ö
get_maximum(xRK4, 0 , 340) = ç	0.064517539485313	÷
ç		÷
è	30	ø

	H (t) := 2	K := 4 m := 8		f(t) := 10
	X0 := 0	V0 := 0	t0 := 0	T := 34 Шаг сетки: h := 0.1
	Решение задачи IV:
	Число узлов сетки:		N :=	T - t0
				h

Формирование вектора правой части системы ОДУ (1) и вектора начальных условий для применения встроенной функции rkfixed:

f1(t , x1 , x2) := x2

f2(t , x1 , x2) :=

f(t) - H(t) ×x2 - K×x1

m

æX0 ö

æf1(t , x0 , x1) ö

x :=

çV0

÷

D (t , x) := çf2

(t , x , x ) ÷

è

ø

è

0 1

ø

10

f(z)

xRK4 := rkfixed(x, t0 , T , N , D)

График решения:

i := 0 .. N

4

xRK4i , 1 2

0 0

20

xRK4i , 0

	9.99
9.98
		0	20		40		60
					z
							æ		0 ö
get_minimum (xRK4, 0 , 340) = ç									0	÷
							ç			÷
							è		0	ø
							æ			4.5	ö
get_maximum(xRK4, 0 , 340) = ç3.92190222525123											÷
							ç				÷
							è			45	ø