66
14.
Заводы							Стоимость перевозки 1т
							1					2				3				4
		Производство це-					цемента, руб.
		мента (т/сут.)					Комбинат					Комбинат				Комбинат					Комбинат
			10				29					50				25					30
1
2			40				20					30				30					22
			10				15					10				30					70
3
4			15				10					15				25					40
			35				20					30				30					10
5
Потребности в цементе							50					20				30					10
15.
Заводы		Производство				Стоимость перевозки 1т цемента, руб.
						Комбинат				Комбинат					Комбинат				Комбинат					Комбинат				Комбинат
		цемента (т/сут.)
						1				2					3				4					5				6
	40						20					5				25				30					25		40
1
	75						210					310				30				22					13		10
2
	60						15					25				70				270					21		12
3
Потребности в цементе							50					20				30				25					30		20

ЗАДЕНИЕ 11 Метод Рунге-Кутты решения дифференциальных уравнений

Решить дифференциальное уравнение численным методом Рунге-Кутты четвертого порядка и построить график с помощью программы Excel.

Дано:

dy = e-nt + mt - ( n + m ) y , dt

где n - номер варианта,

m - двузначный номер группы,

нулевые начальные условия.

Данный метод является одним из методов повышенной точности. Пусть на от-

резке [0,b] требуется найти решение уравнения	dy	= f ( t , y )	(1)

dt

67

с начальным условием y( 0 ) = y0 . Разобьем отрезок [0,b] на n равных частей точ-

ками ti = i × h, ( i = 0 ,1,...,n ) , где h = b / n -

шаг

интегрирования. Каждое после-

дующее значение функции y определяется через предыдущее

yi +1

= yi + Dyi .

(2)

При этом для нахождения приращения Dyi

вычисляются четыре коэффициента

k1( i ) = h × f ( ti , yi );

( i )

æ

h

( i )

ö

k

= h ×

ç t

+

, y

+

k1

÷

;

2

i

i

ç

2

2

÷

è

ø

æ

h

( i )

ö

(3)

k

( i )

= h ×

ç t

+

, y

+

k2

÷

;

3

i

i

ç

2

2

÷

è

ø

k4( i ) = h × f ( ti + h, yi + k3( i ) ).

Затем определяется Dy по формуле

Dyi =

1

( k1( i )

+ 2k2( i ) + 2k3( i ) + k4( i ) )

(4)

6

и подставляется в выражение (2). Вычисление начинается с y0 и заканчивается yn .

Таким образом на каждом шаге необходимо вычислить четыре значения функции

f ( t , y ) , входящих в коэффициенты уравнений (3) сомножителями к шагу интег-

рирования h . Геометрическая интерпретация коэффициентовk1 ,k2 ,k3 ,k4 сле-

дующая. Пусть кривая М0СМ1 (рис.11.1) есть решение дифференциального урав-

нения (1). Точка С этой кривой лежит на прямой, параллельной оси y и делящей отрезок [ ti ,ti +1 ] пополам, B и G – точки пересечения касательной, проведенной к

кривой в точке М0 , с ординатами АС								иN1M1. Тогда число k1 с точностью до
множителя h есть угловой коэффициент касательной кривой М0СМ1, то есть
		k1		= h	dy		= h × tga1 .
					dt
æ			h			k1		ö
ç	+			; yi +					÷	k2 с точностью до
Точка B имеет координаты ç ti		2				2			÷ . Следовательно
è								ø

множителя h есть угловой коэффициент касательной проведенной в точке В (BFотрезок касательной)

68

k2

= h× tga2

Через точку М0

проведем прямую, параллельную отрезку BF, тогда точка D имеет

æ

+

h

; y

+

k

2

ö

и k

= h × tga

есть угловой коэффициент касатель-

координаты ç t

i

2

i

÷

3

3

è

2

ø

ной проведенной к интегральной кривой в точкеD. DR1 – отрезок этой касатель-

ной.

y

R2

α4

M1

R1

D

α3

F

C

α2

G

M0

α1

B

A

N1

t

.

ti

ti+h/2

ti+h

коэффициентовki

Рис. 11.1 Геометрическая

интерпретация

Через

точкуM0

проведем прямую, параллельную DR1 , которая пересечет вер-

тикаль в конце шага в точке R ( t

i

+ h ; y

i

+ k

3

) . Тогда k

4

с точностью до множи-

2

2

теля h определяет угловой коэффициент касательной, проведенной к интеграль-

ной кривой в точке R2 , k4 = h ×tga4 .

Таким

образом, расчетное

значение Dyi

связано с углами a j

j = 1,2 ,3,4 соот-

ношением:

Dyi = 1 h( tga1( i ) + 2tga2( i )

+ 2tga3( i ) + tga4( i ) ).

6

Предварительная

оценка

погрешности

делается после

двойного просчета

по

формуле:

yi* - y( ti

) »

y*

- y

;

i

i

15

69

где y( ti ) - значение точного решения, y*i ; yi - приближенные значения, получен-

ные с шагомh и h / 2 соответственно. Если e - заданная точность, то шаг выби-

рают так, чтобы выполнялось условие: h4 < e .

Метод может быть применен к решению уравнений более высокого порядка -пу

тем сведения к системе дифференциальных уравнений первого порядка. Также

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

Пример решения уравнения

dy = e -t + 120 × t dt

с помощью метода Рунге-Кутта:

program Rung_Kutt; const

N = 100; type

Mt=array[0..N] of real; var

F1:text;

F2:text;

I:integer;

Dy,Dt,K1,K2,K3,K4:real; T,Y,Yst: mt;

function Dy_dt(Tay,Yt:real) : real; begin

Dy_dt:=exp(-Tay)+120.*Tay; end;

procedure Rk; begin

Yst[I]:=Dy_dt(T[I],Y[I]);{ Необязательное определение} K1:=Dt*Dy_dt(T[I],Y[I]); K2:=Dt*Dy_dt(T[I]+Dt/2,Y[I]+K1/2); K3:=Dt*Dy_dt(T[I]+Dt/2,Y[I]+K2/2); K4:=Dt*Dy_dt(T[I]+Dt,Y[I]+K3); Dy:=1/6.*(K1+2*K2+2*K3+K4);

end;

begin Assign(F1,'r.dat'); Assign(F2,'r.rez'); Reset(F1); Rewrite(F2);

		70
Readln(F1,T[0],Y[0],Dt);
for I:=1 to N do
begin
T[I]:=T[I-1]+Dt;
Y[I]:=0.;
Yst[I]:=0.;
end;
for I:=0 to N-1 do
begin
Rk;
Y[I+1]:=Y[I]+Dy;
end;
Yst[N]:=Dy_dt(T[N],Y[I]);
Writeln(F2,' I T	Y	Yst');
for I:=0 to N do
Writeln(F2,I:4,T[I]:8:4,'	',Y[I]:8:4,' ',Yst[I]:8:4,' ');
Close(F1);

Close(F2);

end.

	140
	120
	100
Y,Yst	80						Y
	60						Yst
	40
	20
	0
	0	0,2	0,4	0,6	0,8	1	1,2
				T
.
При открытии файла R.rez с помощью программы Excel необходимо указать в ка-
честве разделителя целой и дробной части точку"." вместо запятой ",".

71