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

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

Федеральное государственное автономное

образовательное учреждение

высшего профессионального образования

«СИБИРСКИЙ ФЕДЕРАЛЬНЫЙ УНИВЕРСИТЕТ»

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

институт

Информационные системы

кафедра

Контрольная работа

по Математическому моделированию

наименование дисциплины

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

Красноярск 2014

1. Цель работы: Изучение методов решения обыкновенных дифференциальных уравнений

Задание: Задание 1. Для охлаждения микропроцессора используется металлический теплоотводящий радиатор. Процесс передачи тепла от радиатора в окружающий воздух описывается дифференциальным уравнением

,

где m и c - масса и удельная теплоемкость материала радиатора, T - температура радиатора, t- время, P - выделяемая микропроцессором мощность,

- отводимое тепло, - коэффициент теплоотдачи конвекцией,
S - площадь поверхности радиатора, Т_С - температура окружающей среды.

Радиатор снабжен вентилятором, который автоматически включается если температура процессора и радиатора превышает допустимый предел, то есть Т>Т_max, и останавливается, если Т<Т_min. Включение обдува эквивалентно изменению коэффициента теплоотдачи по следующему закону:

где 0 - коэффициент теплоотдачи при выключенном вентиляторе, 1 - коэффициент теплоотдачи при обдуве.

Рассчитайте участок зависимости Т(t), на котором система охлаждения выходит на рабочий режим T_min<Т(t)<T_max. Параметры радиатора:

c = 950 Дж/кг·К, m = 0,05 кг, S = 0,04 м². Начальную температуру процессора примите равной Т(t=0) = Т_c = 293 K. Прочие данные указаны в таблице.

Для решения данной воспользуемся методом Эйлера, модифицированным методом Эйлера и методом Рунге-Кутта.

2. Описания методов решения:

2.1 Метод Эйлера

Расчетную формулу метода Эйлера можно получить, используя разложение функции u(x) в ряд Тейлора в окрестности некоторой точки x_i:

. (1)

Если приращение h мало (то есть h<<x_i), то члены ряда, начиная со слагаемого, включающего h во второй степени, могут быть отброшены как малые величины. Тогда из (3)в первом приближении получим

. (2)

Воспользуемся формулой (4), применив ее к единственной известной из условия задачи точкеx₀. Найдем в x₀ производную du(x₀)/dx, подставив (2) в (1):

.

Подставив последнее выражение в (4) и полагая x_i= x₀, получим

или, сокращая обозначения, в окончательном виде

.

Таким образом, (4) при известном значении функции u₀= u(x₀) в начальной точке x₀ позволяет найти приближенное значение u₁= u(x₁) при малом смещении h от x₀. На рис. 1 графически показан начальный шаг решения методом Эйлера.

Рис. 1. Метод Эйлера

Решение можно продолжить, используя найденное значение функции u₁ для вычисления следующего значения - u₂. Распространяя эти рассуждения на последующие точки, запишем расчетную формулу метода Эйлера в виде

. (5)

Из рис. 1 видно, что ошибка метода Эйлера на шаге связана с используемой линейной аппроксимацией u(x). Хотя тангенс угла наклона касательной к кривой точного решения в точке (x₀,u₀) известен и равен du(x₀)/dx, он изменяется при смещении от x₀ до x₁. Следовательно, при сохранении начального наклона касательной на всем интервале h расчет u₁выполняется с погрешностью.

Ошибка метода Эйлера на каждом шаге имеет порядок h², так как члены, содержащие h во второй и более высоких степенях, отбрасываются - см. (3) и (4). Уменьшая h можно снизить локальную ошибку на шаге.

2.2 Модифицированный метод Эйлера

Точность метода Эйлера можно существенно повысить, улучшив аппроксимацию u(x) на рассчитываемом шаге. Для этого при разложении u(x) в ряд Тейлора учтем дополнительно слагаемое, содержащее h² и d²u(xi)/dx² в (3). Определим вторую производную, аппроксимировав ее конечной разностью:

, (6)

гдеДx= h, и .

Подставляя полученное выражение в (3) и отбрасывая члены ряда,
начиная со слагаемого, содержащего h³, запишем

.

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

. (7)

Соотношение (7) дает решение для ui+1 в неявном виде, поскольку ui+1 присутствует одновременно в левой и правой его частях. Следует отметить, что использование неявных методов оправдано тем, что они, как правило, более устойчивы, чем явные. дифференциальное уравнение эйлер рунге

Графически модифицированный метод Эйлера представлен на рис. 2. Из рис. 2 видно, что поправка, учитывающая изменение наклона кривой u(x) заметно уменьшает ошибку на шаге h.

Модифицированный метод Эйлера обеспечивает второй порядок точности. Ошибка на каждом шаге при использовании этого метода пропорциональна h³. Повышение точности достигается за счет дополнительных затрат машинного времени при расчете каждого шага.

Рис. 2. Модифицированный метод Эйлера

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

В модифицированном методе Эйлера для получения второй производной d²u(x_i)/dx² используется конечно-разностная формула (6), включающая значения первой производной u'(x_i) и u'(x_i+h) в начальной и конечной
точках шага. Если подобным же образом вычислить третью производную, рассчитав предварительно вторую производную в двух точках шага, то можно с помощью (3) построить расчетную формулу метода третьего порядка точности. Для этого потребуется определить первую производную u'(x) в дополнительной промежуточной точке между x_i и x_i + h.

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

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

, (8)

где A_n - коэффициенты разложения, k_n - последовательность функций

(9)

б_n,в_nm, 0 <m<n ? p - некоторые параметры.

Неизвестные параметрыA_n,б_n ив_nmможно выбрать из условия

, (10)

где функция ш(h)= u(x_i+h) - о(x_i,h) показывает отклонение приближенного решения о(x_i,h) от точного u(x_i+h). Увеличение параметраp в (8) позволяет сделать погрешность, связанную с заменой точного решения приближенным, как угодно малой.

Предположим, что p = 1. Тогда, подставляя (8) в (10), из условия
ш(0) = ш'(0) = 0 получим A₁ = 1 и ш''(0) ? 0, откуда

,

что соответствует формуле Эйлера (5). Таким же образом можно получить формулы более высоких порядков точности, которые называют методами Рунге-Кутта.

Одним из наиболее известных является вариант метода Рунге-Кутта, соответствующий p = 4. Это метод четвертого порядка точности, для которого ошибка на шаге имеет порядок h⁵. Его расчетные формулы имеют следующий вид:

,

где

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

3. Решение задачи в MathCAD

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

где m иc масса и удельная теплоемкость материала радиатора ,

T температура радиатора ,

t время ,

P выделяемая микропроцессором мощность

Радиатор снабжен вентилятором , который автоматически включается если температура процессора и радиатора превышает допустимый предел, то есть Т >Т max, и останавливается , если Т <Т min. Включение обдува эквивалентно изменению коэффициента теплоотдачи по следующему закону :отводимое тепло , коэффициент теплоотдачи конвекцией , S площадь поверхности радиатора , Т С температура окружающей среды.

Рассчитайте участок зависимости Т (t), на котором система охлаждения выходит на рабочий режим Tmin<Т (t)<Tmax. Параметры радиатора : c = 950 Д ж /к г ·К , m = 0,05 к г , S = 0,04 м 2. Начальную температуру процессора примитеравной Т (t=0) = Т c = 293 K. P=65 Вт , б 0=25 Вт /м3К , б1=82 В т /м 3К , Tmin=313 K, Tmax=353 K.

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