Решение дифференциальных уравнений

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

В Mathcad присутствует много различных функций для решения дифференциальных уравнений. Функция odesolve позволяет записывать уравнение в привычном виде.

Рис. Б.9–Решение дифференциального уравнения с начальными условиями

Функция odesolve требует записи, состоящей из трех частей:

1. ключевое слово given (дано);

2. дифференциальное уравнение и начальные (граничные) условия к нему, граничные условия следует записывать только со штрихом;

3. функция odesolve(x, xk, n), где

x — имя переменной, относительно которой решается уравнение;

xk — конец интервала интегрирования. Начало интервала указано в начальных условиях;

n —параметр, определяющий число шагов интегрирования, на которых вычисляется решение дифференциального уравнения.

Чем больше n, тем больше точность решения уравнения. Параметр n не является обязательным. Если его удалить, Mathcad сам выберет число шагов интегрирования.

Рис. Б.10 – Решение дифференциального уравнения с граничными условиями

Функция odesolve возвращает решение дифференциального уравнения в виде функции, поэтому найденное решение можно интегрировать, дифференцировать (см. рис. б.9 – рис. б.10), и использовать в последующих расчетах.

Функция odesolve решает дифференциальные уравнения как с начальными условиями, так и с граничными. Начальные условия задаются в начале интервала интегрирования, граничные же в двух точках, одна из которых обязательно является началом интервала интегрирования, вторая – произвольная и аргумент её больше, чем в начальной точке. Решение уравнения с начальными условиями показано на рис. б.9, с граничными условиями — на рис. б.10.

Решение системы дифференциальных уравнений показано на рис. б.11. Как видно, обращение к функции odesolve изменилось. Для решения одного уравнения: odesolve(x, xk, n).

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

На рис. б.11 при решении системы трех уравнений вектор имен неизвестных содержит 3 имени.

Рис. Б.11 – Решение системы дифференциальных уравнений с начальными условиями

Б.4 Ввод греческих букв

Для ввода греческих букв в математической панели Mathcad есть панель инструментов Greek, которая содержит все греческие символы (рис. б.2).

Знак бесконечности ∞ можно вставить с математической панели Calculs (рис. б.2).

Часто используемое число π можно набрать с математической панели Calculator, где есть кнопка π, или нажав аналогичную кнопку панели греческих букв (рис. б.2).

Более простой способ ввода греческих букв: надо с клавиатуры набрать латинский аналог греческой буквы (первую букву в латинском названии греческой буквы), затем нажать сочетание клавиш Ctrl+G. На экране появится греческое написание буквы.

Греческие буквы и их аналоги (в скобках):

α(a), β(b), χ(c), δ(d), ε(e), η(h), γ(g), λ(l), μ(m), ν(n), ω(w), ϕ(f), π(p), ψ(y), ρ(r), σ(s), τ(t), θ(q), ξ(x), ς(z).

Если какой-то константе (переменной) не присвоено никакого значения левее и выше ее положения на экране, то она окрашена красным цветом, и это свидетельствует об ошибке.

Б.5 Ввод текста

Для ввода текста в документ Mathcad можно в главном меню выбрать команду Insert Text Region (Вставить Текстовую область), или ввести с клавиатуры символ " (кавычка), при этом на экране появляется текстовая область, в которой можно печатать текст.

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

В текстовую область можно вставлять математическую. В главном меню Mathcad необходимо выбрать команду Insert Math. Region (Вставить Математическую область). Вставленная математическая область участвует в вычислениях наравне с другими выражениями.

При желании вставленную математическую область можно отключить (рис. б.12). Для этого необходимо щелкнуть правой кнопкой мыши на вставленном выражении и в открывшемся меню выберать команду Disable Evaluation (Отключить вычисление). После этого вставленное выражение будет играть роль рисунка.

Рис. Б.12 – Математическая область, вставленная в текст

Текстовой области в Mathcad присвоен стиль Normal. Чтобы настроить его, необходимо в главном меню выбрать команду Style Normal Modify Font (Стиль Обычный Изменить Шрифт) и в появившемся диалоговом окне выбрать необходимые характеристики шрифта.

Шрифты Mathcad плохо воспринимают кириллицу. Привычный по работе в Word шрифт Times New Roman кириллицу не принимает. С кириллицей работают шрифты System и Ms Sans Serif.

ПРИЛОЖЕНИЕ В Варианты заданий

ПРИЛОЖЕНИЕ Г Пример выполнения задания

В качестве примера рассмотрим исходные данные варианта № 30. Схема фильтра, величины сопротивлений, емкостей, индуктивностей и независимые начальные условия для дифференциальных уравнений представлены на рис. г.1.

Рис. Г.1 ‑ Схема, элементы, начальные условия фильтра

Г.1 Частотные характеристики фильтра

Комплексные сопротивления ветвей фильтра ( , , см. табл. а.1 и табл. а.2) вычисляем по следующим формулам

, , .

График передаточной функции фильтра , размещенный на комплексной плоскости, представлен на рис. г.2 АЧХ и ФЧХ передаточной функции представлены на рис. г.3.

Графики, для заданных вариантом задания значений сопротивлений, емкостей и индуктивностей построены средствами пакета MathCad.

Рис. Г.2 ‑ Передаточная функция фильтра

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

Рис. Г.3 – АЧХ и ФЧХ фильтра

Г.2 Система дифференциальных уравнений

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

Г.2.1 Составление системы уравнений

Составим уравнения фильтра по контурам и узлам, показанным на рис. 2.13. Интегро-дифференциальные уравнения, составленные по контурам 1), 2), 3), выглядят следующим образом

,

алгебраическое уравнение, составленное по узлу а), выглядит следующим образом

.

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

.

Зависимые начальные условия определяем из систем уравнений , ,

.

Дифференцируем, левые и правые части составленных интегро-дифференциальных уравнений, и получаем следующую систему линейных неоднородных дифференциальных уравнений с постоянными коэффициентами

.

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

.

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

Г.2.2 Решение системы средствами Odesolve

Для решения полученной системы уравнений используем встроенную процедуру Odesolve.

На рис. г.4 представлена процедура интегрирования полученных уравнений со следующими исходными данными:

- R, C, L, Rn – значения величин элементов;

- t - переменная интегрирования;

- dlt ( ) - шаг интегрирования;

- steps - число шагов интегрирования;

- T1 - время интегрирования;

- e(t) - единичное входное воздействие.

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

Рис. Г.4 – Процедура интегрирования системы Odesolve

Рис. Г.5 – Результаты интегрирования системы. Odesolve

Преобразуем полученную систему уравнений в систему разностных уравнений по методу Эйлера.

Г.2.3 Система разностных уравнений

Заменим производную в нормальной системе дифференциальных уравнениях следующей разностью

,

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

.

Оставим в левой части уравнения старшие (по времени) переменные и соберем в правой части все другие переменные

.

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

Г.2.4 Решение системы разностных уравнений

Процедура интегрирования разностных уравнений Mysolve представлена на рис. г.6. Результат интегрирования уравнений системы Mysolve представлена на рис. г.7.

Рис. Г.6 – Процедура интегрирования системы Mysolve

Рис. Г.7 – Результаты интегрирования системы. Mysolve

Г.2.5 Сравнение полученных решений

Сравним результаты решений системы процедурами Odesolve и Mysolve путем наложения графиков решений, полученных этими процедурами. Отобразим графики переменной , на одной временной оси (см.рис. г.8).

Рис. Г.8 – Сравнение результатов интегрирования системы

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

Г.3 Дифференциальное уравнение 3-го порядка

Г.3.1 Получение дифференциального уравнения

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

.

В этой системе, помимо переменной и ее трех производных, имеются следующие другие переменные

, , , .

Для получения нулевых коэффициентов при других переменных дополним заданную систему уравнений, двумя другими линейно независимыми уравнениями.

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

Четвертое уравнение в новой системе получим, дифференцируя алгебраическое уравнение 3).

Пятое уравнение получим, дифференцируя уравнение 2).

Запишем новую систему уравнений, оставляя в правой части уравнений ноль

.

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

, , ,

, ,

, .

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

,

где коэффициенты , , , , будут следующие

,

,

,

,

.

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

.

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

Г.3.2 Сравнение частотных характеристик

Получим передаточную функцию фильтра , используя дифференциальное уравнение 3-го порядка. Определим нули и полюсы передаточной функции, используя встроенную процедуру Solve. Процедура определения полюсов и нулей и получения передаточной функции представлено на рис. г.9. График полученной функция представлена на рис. г.10. Сравнение частотных характеристик представлено на рис. г.11.

Рис. Г.9 – Получение передаточной функции. Solve

Рис. Г.10 – Передаточная функция уравнения 3-го порядка

Рис. Г.11 – Сравнение АЧХ и ФЧХ фильтра

Частотные характеристики функции, полученной ранее и частотные характеристики и функции , совпадают.

Г.3.3 Решение уравнения средствами Odesolve

На рис. г.12 представлена процедура интегрирования полученных уравнений со следующими исходными данными:

R, C, L, Rn – значения величин элементов;

- t - переменная интегрирования;

- dlt ( ) - шаг интегрирования;

- steps - число шагов интегрирования;

- T1 - время интегрирования;

- e(t) - единичное входное воздействие.

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

Рис. Г.12 – Процедуры интегрирования уравнения Odesolve

Рис. Г.13 – Результаты интегрирования уравнения. Odesolve

Преобразуем, по методу Эйлера, полученное уравнение 3-го порядка в разностное уравнение.

Г.3.4 Разностное уравнение

Заменим производные в дифференциальном уравнении 3-го порядка следующими разностями

,

,

.

Проведя замену, получим следующее разностное уравнение

.

Сгруппируем коэффициенты и получим следующие уравнение

.

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

,

,

,

,

.

Выражение, определяющее зависимость значения переменной на шаге n ( ) от значений этой переменной на предыдущих шагах ( , , ) и значений правой части уравнения ( , ), выглядит следующим образом

.

Используя это уравнение, напишем процедуру Mysolve.

Г.3.5 Решение разностного уравнения

Процедура интегрирования разностного уравнения Mysolve представлена на рис. г.16. Результат интегрирования уравнения процедурой Mysolve представлена на рис. г.14.

Г.3.6 Сравнение полученных решений

Сравним результаты решений уравнения процедурами Odesolve и Mysolve путем наложения графиков решений, полученных этими процедурами. Отобразим графики переменных и , на одной временной оси (см. рис. г.15).

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

Рис. Г.14 – Результаты интегрирования уравнения. Mysolve

Рис. Г.15 – Сравнение результатов интегрирования уравнения

Рис. Г.16 – Процедура интегрирования уравнения Mysolve

Г.4 Программирование в среде Code Vision AVR

Г.4.1 Решение системы по разностной схеме

Программа соответствующая системе разностных уравнений, написанная на языке С++, представлена ниже.

#include <mega8535.h>

#define ADC_VREF_TYPE 0x60

unsigned int ss, dat, YY, Y2;

double Y1; EE_n, EE_n_1, EE_n_2, X1_n, X1_n_1, X1_n_2,

X2_n, X2_n_1, X2_n_2, X3_n, X3_n_1, X3_n_2,

R1=10.0, R2=0.0, Rn=0.9,

C1=0.01, C2=0.08, L3=0.14, ZZ=0,

cnt=0, dt=0.0215;

unsigned char adc_input=0, a;

// задание временных интервалов таймер T0

interrupt [TIM0_OVF] void timer0_ovf_isr(void)

{ TCNT0=235; // начальная инициализация таймера

PORTB.1=!PORTB.1;

ADMUX=adc_input|ADC_VREF_TYPE;

ADCSRA|=0x40; // запуск АЦП }

interrupt [ADC_INT] void adc_isr(void)

{unsigned char adc_data;

adc_data=ADCH; // получение данных напряжения от АЦП

Y1=ADCH;

Y1=Y1/255-0.5*0;

// задание начальных условий

// начальные условия в момент времени 0

if (cnt==0){

EE_n_2=Y1;

X1_n_2=Y1/R1;

X2_n_2=-Y1/R1;

X3_n_2=0;

X1_n_1=0;X2_n_1=0;X3_n_1=0;

X1_n=0;X2_n=0;X3_n=0;

// начальные условия в момент времени 1

PORTA.1=1;}

if (cnt==1){

EE_n_1=Y1;

X1_n_1=X1_n_2+(1/R1)*((EE_n_1-EE_n_2)/dt)*dt

-(1/(R1*C1))*X1_n_2*dt-(1/(R1*C2))*X2_n_2*dt;

X3_n_1=X3_n_2;

X2_n_1=-X1_n_1-X3_n_1;

X1_n=0;X2_n=0;X3_n=0;

PORTA.1=0;}

// решение диф. Уравнения методом Эйлера

if (cnt>1){

EE_n=Y1;

X1_n=X1_n_1+(1/R1)*((EE_n-EE_n_1)/dt)*dt-(1/(R1*C1))*X1_n_1*dt

-(1/(R1*C2))*X2_n_1*dt;

X3_n = 2*X3_n_1-X3_n_2 - (Rn/L3)*((X3_n_1-X3_n_2)/dt)*dt*dt

+ (1/(C2*L3))*X2_n_2*dt*dt;

X2_n=-X1_n-X3_n;

// рекурсия

EE_n_2=EE_n_1; EE_n_1=EE_n;

X3_n_2=X3_n_1; X3_n_1=X3_n;

X2_n_2=X2_n_1; X2_n_1=X2_n;

X1_n_2=X1_n_1; X1_n_1=X1_n;

}

// вывод результатов на виртуальный ЦАП VMLAB

PORTC=X3_n*80*20+255/2; // ток i3

PORTD=X2_n*80*20+255/2; // ток i2

PORTB=X1_n*80*20+255/2; // ток i1

// контроль времени интегрирования

cnt=cnt+1;

if(cnt>100){cnt=0;}

}

void main(void)

{

// инициализация микроконтроллера

// конфигурирование порта A на ввод

PORTA=0x00;

DDRA=0xfe;

// конфигурирование порта B на вывод

PORTB=0x00;

DDRB=0x00;

// конфигурирование порта C на вывод

PORTC=0x00;

DDRC=0xff;

// конфигурирование порта D на вывод

PORTD=0x00;

DDRD=0x00;

// инициализация таймера временных интервалов

TCCR0=0x05; // установка режима таймера

TCNT0=0x00; // установка начальных значений

OCR0=0x00;

// разрешение прерываний о т таймера

TIMSK=0x01; // установка режима таймера T0

// инициализация АЦП

// установка разрядности АЦП

// установка опорного напряжения АЦП

// установка способа запуска АЦП

ADMUX=ADC_VREF_TYPE;

ADCSRA=0x8D;

SFIOR&=0xEF;

// глобальное разрешение прерываний

#asm("sei")

while (1)

{// холостой цикл };

}

Г.4.2 Результаты решения системы

Графики решений разностных уравнений, построенные средствами среды Visual Micro Lab, представлены на рис. г.17.

Рис. Г.17 – Результаты интегрирования. С++

Г.5 Выводы по проделанной работе

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

Полученные частотные характеристики заданного фильтра указывают на то, что это фильтр нижних частот.

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

ания. С++. 162