Informatika_met_ukazania

	Если		обыкновенное дифференциальное		уравнение	2-го порядка
y		d 2 y(x)		f (x, y(x), y '(x)), (x, y(x)) G	R3 имеет	в области G
		dx2

решение, то, вообще говоря, таких решений бесконечно много и они могут быть заданы в виде функций y y(x, C) , где С – произвольная

константа, такая, что (x, y(x, C)) G и y ''(x, C) f (x, y(x, C), y (x,C)) при произвольных значениях С. Если же необ-

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

ловию y(x0 )		y0 , y (x0 ) y10 , то говорят об отыскании решения задачи
Коши.
Численное решение нормальной системы дифференциальных
уравнений 1-го порядка с системой начальных условий
	y	f (x, y , y			2	, ...,	y	n	),	y1 (x0 )	y0,1 ,
	1	1	1
	y2	f2 (x, y1		, y2 , ...,			yn ),			y2 (x0 )	y0,2 ,
	. . . . . . . . . . . . . . . . . . . . . .									. . . . . . . . . . . .
	yn	fn (x, y1		, y2 , ...,			yn ),			yn (x0 ) y0,n
состоит	в	построении					таблицы			приближенных значений
yi,1, yi,2 ,	. . . ,	yi,n ,	i = 1, 2, ..., N, решения							y1 (x), y2 (x), ..., yn (x) на от-

резке [x0, xN] в точках x0, x1, ..., xN, которые называются узлами сетки.

Всистеме MatLab решение дифференциальных уравнений сводится

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

Пример 1. Найдем на отрезке [0, 3] приближенное решение задачи Коши для системы дифференциальных уравнений

y1	y2 sin(xy3 ),	y (0) 1,
		1

y	2	y2	,	с начальными условиями	y	2	(0)	0,
		1
y3		y3 y1			y3		(0)	1

и построим графики для найденного решения.

41

На первом этапе построим специальную М-файл функцию для вычисления правых частей системы. Она должна иметь два входных аргумента: независимую переменную, по которой производится дифференцирование (в нашем случае, очевидно, ею является х − даже если х явно не входит в систему), и вектор искомых функций (в нашем случае, очевидно, векторстолбец y с тремя элементами). Выходным аргументом М-файл функции является вектор правой части системы − как показано на рис. 4.

Рис. 4. Окно редактора М-файлов с функцией dif3

На втором этапе вызывается один из семи солверов (решателей) системы MatLab − мы ограничимся использованием лишь одного из них − ode45, который реализует одношаговый метод Рунге–Кутта четвѐртого порядка. Он рекомендуется разработчиками MatLab в качестве первого для использования. Решение находится в виде таблицы чис-

ленных значений искомых функций y1 (x), y2 (x), ..., yn (x) для некото-

рого набора значений аргумента x.

До вызова солвера необходимо сформировать вектор-столбец начальных условий, например:

>> Y0=[1; 0; 1];

Синтаксис вызова солвера − один из трех следующих:

[X,Y]=ode45(‘fun’,[x0 xE],y0), [X,Y]=ode45(‘fun’,[x0:step:xE],y0), [X,Y]=ode45(‘fun’,[x0,x1,x2,…,xE],y0),

где Х − вектор-столбец значений независимой переменной с начальным значением х0 и конечным значением хЕ; step − шаг изменения этой независимой переменной; x0, x1, x2, …, xE − конкретная возрастающая или убывающая числовая последовательность, для которой необходимо найти искомое решение; fun − имя М-файла (в нашем случае dif3); y0 − имя вектора начальных условий.

42

Пример 1.

>> [T,Y]=ode45('dif3',[0:0.5:3],Y0) T =

0
0.5000
1.0000
1.5000
2.0000
2.5000
3.0000
Y =
1.0000	0	1.0000
0.9324	0.4826	0.2211
0.5914	0.7973	-0.1700
0.0321	0.8621	-0.2202
-0.5190	0.9069	-0.0249
-0.8896	1.1735	0.2751
-1.1528	1.6915	0.5705

Вычисленные значения искомых функций (в нашем случае это y1 (x) , y2 (x) и y3 (x)) располагаются по столбцам.

На последнем, третьем, этапе по значениям вектора-столбца Т и столбцам матрицы Y можно, очевидно, построить графики (рис. 5).

>>plot(T,Y(:,1),'k-o',T,Y(:,2),'b-s',T,Y(:,3),'r-v')

>>grid on

Рис. 5. Окно выходного графика

43

Пример 2. Найдем на отрезке [0, ] приближенное решение ОДУ 1-го порядка y sin(xy), удовлетворяющее начальным условиям

y(0) = 1, и построим график найденного решения. Шаг изменения независимой переменной х − 0.4.

Вид М-файла с функцией dif1 представлен на рис. 6, график искомой функции − на рис. 7, а в командном окне необходимо, например, подать следующие команды:

>>Y0=[1];

>>[T,Y]=ode45('dif1',[0:0.4:pi],Y0);

>>plot(T,Y(:,1),'k-o')

>>grid on

Рис. 6. Окно редактора М-файлов с функцией dif1

Рис. 7. Окно выходного графика

44

В радиотехнике ОДУ 1-го порядка применяются для анализа простейших электронных цепей во временной области, например, показанных на рис. 8. Там же приведены аналитические выражения, описывающие их. Искомыми при решении задачи анализа могут быть следующие параметры: i(t), uR (t) Ri(t) , uL (t) , uC (t) . Очевидно, что

для решения задачи Коши необходимы также и начальные условия − для цепи, приведенной на рис 8,а, их обычно задают в виде i(t0 ) i0 ,

а для цепи (рис. 8,б) − uC (t0 ) uC0 .

а

б

Рис. 8. Простейшие электронные цепи и описывающие их уравнения

Для составления дифференциальных уравнений, описывающих электронные цепи, используют законы Кирхгофа, известные из курса физики. При этом напряжения и токи на резисторе R, индуктивности L и емкости C связаны соответственно следующими выражениями:

u (t) Ri(t) , u

(t) L

di(t)

и

u (t)

1

i(t)dt ,

u (t)

1

i(t) .

L

R

dt

C

C

C

C

Пример 3. Для цепи (рис. 8,а)

построить

зависимости

uR (t) и

uL (t)

для диапазона времен 0

t 120 мкс при R 5 Ом, L

0.1 мГн

и ЭДС в виде ступеньки напряжения

e(t)

0 при t

5 мкс,

1 при t

5 мкс.

45

Число узлов сетки выбрать равным 100. Очевидно, что при t0 0 мы имеем i(t0 ) 0 − это и будут наши начальные условия. Фрагмент

соответствующего командного окна, решающего поставленную задачу, а также виды окна М-файла и итогового графического окна представлены ниже и на рис. 9 и 10.

Y0=[0];	% начальные условия
[T,Y]=ode45('rl',[0:120e-6/100:120e-6],Y0);
UR=5*Y;	% напряжение на резисторе
for i=1:1:101

T1(i)=120e-6/100*(i-1);

if T1(i)<=5e-6
E(i)=0;
else
E(i)=1;
end
end
UL=E'-UR;	% напряжение на индуктивности

plot(T,Y(:,1),'k--',T,UR(:,1),'r-.',T1,UL(:,1),'b-',T1,E,'m-') grid on

8

Рис. 9. Окно редактора М-файлов с функцией rl

Для того чтобы в системе MatLab решить задачу Коши для дифференциального уравнения n-го порядка

	y n	d n y(x)	f (x, y(x), y '(x), ..., y n 1	(x)), (x, y(x)) G R3
		dxn
			46

с начальными условиями y(x )	y , y (x )		y	, ..., y(n 1) (x ) y	n 10	,
0	0	0	10	0

его приводят, как уже было сказано, к системе из n уравнений первого порядка и далее следуют описанной ранее схеме из трех этапов.

Рис. 10. Окно выходного графика

Пример 4. Пусть требуется решить дифференциальное уравнение

2-го	порядка	вида y 2y 20y cos(x) при начальных условиях
y(0)	1, y (0)	0 на отрезке [0, 10] с шагом 0.1. Для приведения диф-

ференциального уравнения n-го порядка к системе из n дифференциальных уравнений 1-го порядка вводят вспомогательные функции. Для

	нашего случая это выглядит так: y1 y , y2			y . После этого диффе-
	ренциальное уравнение можно заменить системой
	y1	y2 ,	при начальных условиях		y1 (0)	0,
	y2	2 y2 20 y1 cos(x)			y2 (0)	1.

Для решения этой системы дифференциальных уравнений следует создать М-файл функцию, показанную на рис. 11, и в командном окне набрать следующий код (результат виден из рис. 12):

>>Y0=[0;1];

>>[T,Y]=ode45('dif2',[0:0.1:10],Y0);

>>plot(T,Y(:,1),'k--',T,Y(:,2),'r-')

>>grid on

47

Рис. 11. Окно редактора М-файла с функцией dif2

Первым столбцом найденной матрицы Y являются значения искомой функции y(x) , вторым столбцом − значения ее производной y (x) (обратите внимание на рис. 12).

Рис. 12. Окно выходных графиков

В радиотехнике ОДУ 2-го порядка применяются для анализа простых электронных цепей во временной области, например, показанных на рис. 13. Там же приведены аналитические выражения, описывающие их. Искомыми при решении задачи анализа могут быть следующие параметры: i(t), uR (t) Ri(t) , uL (t) , uC (t) (для рис. 13,а) и

u(t), iR (t) u(t) / R , iL (t) , iC (t) (для рис. 13,б). Очевидно, что для решения задачи Коши необходимы также и начальные условия − для це-

пи, представленной на рис. 13,а, их обычно задают в виде i(t0 ) i0 , i (t0 ) i10 , а для цепи (рис. 13,б) − u(t0 ) u0 , u (t0 ) u10 .

48

а

б

Рис. 13. Простые электронные цепи и описывающие их уравнения

Пример 5. Для цепи (рис. 13,а) построить зависимости uR (t) , uL (t) и uC (t) в диапазоне времен 0 t 25 мкс при R 5 Ом, L 10 мкГн ,

С2.533 нФ и ЭДС в виде ступеньки гармонического колебания

e(t)	0 при t 5 мкс,	,
	cos(2 ft) при t 5 мкс,

49

где f 1 МГц . Шаг сетки времени выбрать равным 10 нс. Очевидно,

что при t0 0 мы имеем i(t0 ) 0 , i (t0 ) 0 − это и будут наши начальные условия. Фрагмент соответствующего командного окна, ре-

шающего поставленную задачу, а также вид окон М-файлов и итогового графического окна (только для напряжения на резисторе) представлен ниже и на рис. 14 и 15.

>>Y0=[0;0];

>>[T,Y]=ode45('rlcpos',[0:0.1e-7:25e-6],Y0);

>>plot(T,5*Y(:,1)) % напряжение на резисторе

>>grid on

а

б

Рис. 14. Окна редакторов М-файлов с функциями

50