Применение системы MathCAD для исследования температуры тела

  Имеется два типа задач, для которых возможно численное решение ОДУ с помощью MathCAD:

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

Из дифференциальных уравнений в частных производных есть возможность решать только уравнения с двумя независимыми переменными: одномерные параболические и гиперболические уравнения, такие как уравнения теплопроводности, диффузии, волновые уравнения, а также двухмерные эллиптические уравнения (уравнения Пуассона и Лапласа).

В MathCAD нет универсальной функции для решения дифференциальных уравнений, а есть около двадцати функций для различных видов уравнений, дополнительных условий и методов решения. Эти функции можно найти в библиотеке Insert/Function, категория “Differential Equation Solving (решение дифференциальных уравнений).

 

Решение Обыкновенных Дифференциальных Уравнений (ОДУ)

ОДУ первого порядка.

ОДУ первого порядка  называется уравнение

F(x,y,y’)=0

где: F – известная функция трех переменных;

x – независимая переменная на интервале интегрирования[a,b];

y – неизвестная функция;

y’ – ее производная.

Функция y(x) является решением дифференциального уравнения, если она при всех xÎ[a,b] удовлетворяет уравнению

F(x,y(x),y’(x))=0

График решения y(x) называется интегральной кривой дифференциального уравнения. Если не заданы начальные условия, таких решений y(x) будет множество. При известных начальных условиях y(x₀)= y₀ решение y(x) будет единственным.

Вычислительный процессор MathCAD может работать только с нормальной формой ОДУ.

  Нормальная форма ОДУ – это ОДУ, разрешенное относительно производной

y’=f(x,y)

ОДУ высших порядков

Обыкновенным дифференциальным уравнением n-го порядка называется уравнение вида:

F(x,y,y’,y’’, …,y⁽ⁿ⁾)=0

где:    F – известная функция n+2 переменных;

x – независимая переменная на интервале интегрирования[a,b];

y – неизвестная функция;

n – порядок уравнения.

Функция y(x) является решением дифференциального уравнения, если она при всех xÎ[a,b] удовлетворяет уравнению

F(x, y(x), y’(x), y’’(x),…, y⁽ⁿ⁾(x))=0

Нормальная форма ОДУ высшего порядка имеет вид:

Y⁽ⁿ⁾ =f(x, y, y’, …, y^(n-1))

Если не заданы начальные  условия, то дифференциальное уравнение n – го порядка имеет бесконечное множество решений, при задании начальных условий y(x₀)= y₀, y’(x₀)= y_0,1, y’’(x₀)= y_0,2, …, y^(n-1)(x₀)= y_0,n-1 решение становится единственным (задача Коши).

Задача Коши для дифференциального уравнения n – го порядка может быть сведена к задаче Коши для нормальной системы n дифференциальных уравнений 1-го порядка, которая в векторной форме имеет вид

Y’ = F(x, Y),     Y(x₀) = Y₀

где:    Y(x₀) = Y₀ – вектор начальных условий;

Y’=(y’₁, y’₂, …, y’_n) – вектор первых производных;

F(x, Y) = (y₂, y₃, …, y_n, f(x,y₁, … , y_n) – вектор правых частей;

Y = (y₂, y₃, …, y_n) – вектор искомого решения.

Эта система получается в результате следующей замены: 

 

 

 

 

 

 

 

Для численного интегрирования ОДУ в MathCAD имеется выбор – либо использовать вычислительный блок Given/Odesolve, либо встроенные функции. Оба способа обладают одинаковыми возможностями, но при использовании блока решения запись уравнений более привычна и наглядна, однако отдельная функция может быть использована в составе других функций и программ. Рассмотрим оба варианта решения.

Для численного интегрирования ОДУ в MathCAD имеется выбор – либо использовать вычислительный блок Given/Odesolve, либо встроенные функции. Оба способа обладают одинаковыми возможностями, но при использовании блока решения запись уравнений более привычна и наглядна, однако отдельная функция может быть использована в составе других функций и программ. Рассмотрим оба варианта решения.

Вычислительный блок Given/Odesolve    

Ниже приведены два  примера для решения дифференциальных уравнений первого и второго  порядка с использованием вычислительного  блока решения Given/Odesolve

 

1.3 РЕАЛИЗАЦИЯ ЧИСЛЕННЫХ МЕТОДОВ

Вычислительный блок для решения одного ОДУ состоит  из трех частей:

• ключевое слово given;
• ОДУ и начальные условия, записанные с помощью логического равенства;
• встроенная функция Odesolve(x, b) относительно независимой переменной x на интервале [a, b]; b – верхняя граница отрезка интегрирования.

Допустимо и даже предпочтительнее задание функции Odesolve(a, b, step) с тремя параметрами, где step – внутренний параметр численного метода,

определяющий количество шагов; чем больше step, тем с лучшей точностью

будет получен результат, но тем больше времени будет затрачено на его поиск.

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

Функция Odesolve использует для решения дифференциальных уравнений наиболее популярный алгоритм Рунге-Кутта четвертого порядка, описанный в большинстве книг по методам вычислений. Он обеспечивает малую погрешность для широкого класса систем ОДУ за исключением жестких систем. Если щелчком правой кнопки мыши на блоке формул с функцией Odesolve вызвать контекстное меню, то можно изменить метод вычисления решения, выбрав один из трех вариантов: Fixed – метод Рунге-Кутта с фиксированным шагом интегрирования (этот метод используется по умолчанию), Adaptive – также метод Рунге-Кутта, но с переменным шагом, изменяемым в зависимости от скорости изменения функции решения, Stiff – метод, адаптированный для решения жестких уравнений и систем (используется так называемый метод PADAUS).

Альтернативный метод решения ОДУ заключается в использовании одной из встроенных функций: rkfixed, Rkadapt, или Bulstoer. Все они решают задачу Коши для системы дифференциальных уравнений первого порядка, но каждая из них использует для этого свой метод. Для простых систем не играет большой роли, какой метод использовать – все равно получите решение достаточно быстро и с высокой точностью. Но для сложных или специфических систем бывает, что некоторые методы вообще не могут дать удовлетворительного решения за приемлемое время. Именно для таких сложных, но не редких случаев в MathCAD и введено несколько различных методов решения систем ДУ.

• rkfixed – метод Рунге-Кутта с фиксированным шагом интегрирования. Самый простой и быстрый метод, но далеко не всегда самый точный. Полностью аналогичен использованию функции Odesolve с выбранным в контекстном меню методом Fixed.
• Bulstoer – метод Булирша – Штера. Этот метод более эффективен, чем метод Рунге-Кутта, в случае если решение является плавной функцией.

Имена функций Rkadapt и Bulstoer начинаются с прописной буквы. В MathCAD для некоторых имен функций неважно, с какой буквы они записаны, но для перечисленных функций это принципиально, т.к. в MathCAD также существуют функции с такими же именами, только записанные с маленькой буквы – rkadap, bulstoer. Эти функции используются в тех случаях, когда важным является решение задачи в конечной точке интервала интегрирования.  

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

Рисунок 1.4 – Решение ОДУ первого порядка

 

 

Рисунок 1.5 – Решение ОДУ высшего порядка

 

 

 

 

 

2 АЛГОРИТМИЧЕСКИЙ  АНАЛИЗ В ЗАДАЧИ

1. ПОСТАНОВКА ЗАДАЧИ

1. Тело, имеющее в начальный момент температуру T(0)=T₀, поместили в среду, температура которой поддерживается постоянной и равна T.

С использованием системы MathCAD определить изменения температуры тела с течением времени. Построить график функции T(t), сделать выводы.   К чему стремится значение функции T(t) с увеличением t?

2. Исследовать влияние значений изменяемого параметра на вид функции температуры в среде. Значения варьируемого параметра выбрать самостоятельно (не менее 7 значений).

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

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.2 АНАЛИЗ ИСХОДНЫХ  ДАННЫХ. ОПИСАНИЕ МАТЕМАТИЧЕСКОЙ МОДЕЛИ

Исходными данными для  работы являются:

T₀ – начальная температура тела;

T₁ –  температура среды; 

t – время исследования;

 – коэффициент пропорциональности ( > 0);

Таблица 2.1 – исходные данныx:

Вариант	T0, С0	T1, С0		t, с	Варьируемый параметр
3-6	50	30	0.087	0-150	T1

 

 

 

 

 

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

Тело, имеющее в начальный момент температуру T(0)=T₀, поместили в среду, температура которой поддерживается постоянной и равна T.

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

где коэффициент пропорциональности.

Знак минус в правой части  уравнения соответствуют экспериментальным  данным: если T(t)-T₁>0, то температура тела убывает и поэтому скорость его изменения отрицательна, если же T(t)-T₁<0, то температура тела возрастает, а следовательно, скорость ее изменения положительна.

Итак, процесс нагревания (или охлаждения) тела в среде с изменной температурой моделируется уравнением указанным выше.

 

 

Значения варьируемого параметра T₁ выбираем самостоятельно. Значения варьируемого параметра запишем в таблицу 2.

Таблица 2.2 – Значения варьируемого параметра T₁

T1 С0	T2 С0	T3 С0	T4 С0	T5 С0	T6 С0	T7 С0	T8С0
0	5	10	15	20	25	30	35

 

Описание математической модели

1. Задаём начальные условия:

T₁=32;  y=0.087.

2. Записываем исходное уравнение:

D (t,T):=-y(T-T₁)

3. Решаем уравнение с помощью функции rkfixed, получим:

Z:=rkfixed (T₀,0,65,1000,D)

4. Получаем значения Z, где 1-ый столбец – это время, а второй – температура в определённое время.

Рисунок 2.1 – Таблица значений

5)  После этого, строим график зависимости температуры от времени, а так же определим изменение температуры в течении определённого времени.

Рисунок 2.2 – График зависимости

6)  Для определения изменения температуры в течении определённого времени: в таблице Z находим значение t, равное температуре T.

Изменяем значение параметра не менее 7 раз.

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

2.3 СХЕМА АЛГОРИТМА  РЕШЕНИЯ ЗАДАЧИ И ЕЁ ОПИСАНИЕ

 

 

 

 

 

 

 

 

 

1. В MathCAD вводим исходные данные.
2. Нашли  функции зависимости температуры тела от времени.
3. Построили  график зависимости температуры тела от времени.
4. Построили сводный график зависимостей температуры тела от времени.
5. Подобрали  аппроксимирующую функцию, по результатом исследования.
6. Построили график исходной и аппроксимирующей функций .

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3.ОПИСАНИЕ  РЕАЛИЗАЦИИ ЗАДАЧИ

3.1ОПИСАНИЕ РЕАЛИЗАЦИИ  ЗАДАЧИ В MathCAD

Задаём начальные условия. Известно, что скорость изменения  температуры тела пропорциональна  разности температур тела и окружающей среды и выражается уравнением: