3.2. Метод эйлера-коши
Пусть опять решаем уравнение y'= f(x,y),
y(x0)= y0.
Решение ищем на отрезке [x0,xn].
Пусть нам известны координаты некоторой точки, принадлежащей искомому решению (xm,ym). Найдём средний тангенс угла наклона касательной для двух точек: (xm,ym) и (xm+h,ym+h·ỹm).
Последняя точка, есть та самая, которую в методе Эйлера мы обозначаем (xm+1,ym+1), но здесь точка будет вспомогательной. ỹm+1
ỹm+1 ym+1 ym
Итак, сначала по методу Эйлера находится точка А, лежащая на прямой L1, тангенс угла наклона которой tgα1= f(xm,ym). В этой точке снова вычисляется тангенс угла наклона касательной L2
tgα2= f(xm+h,ym+hỹm)
Затем через точку (xm,ym) проводим прямую L, тангенс угла наклона которой равен (tgα1+tgα2)/2. Точка, в которой L пересекается с прямой x=xm+1 и будет искомой (xm+1,ym+1). Таким образом, ym+1 есть искомое приближение значений функции на данном шаге интегрирования.
Расчетные формулы метода Эйлера-Коши следующие:
ỹm+1= ym+hf(xm,ym);
ym+1= ym+(h/2)·[f(xm,ym)+f(xm+1,ỹm+1)];
xm+1= xm+h.
Аналогично, для системы дифференциальных уравнений:
yi'= fi(x,y1,y2,…yk), yi(x0)=yio, i=1,2,..k;
ỹim+1= yim+ hfi(xm,y1m,y2m,…ykm);
yim+1= yim+ (h/2)·[fi(xm,y1m,y2m,…ykm)+fi(xm+1,ỹ1m+1,ỹ2m+1,ỹkm+1)];
xm+1= xm+h.
Здесь i – номер уравнения системы, m – номер шага.
ПРОГРАММА МЕТОДА ЭЙЛЕРА.
PROGRAM EILER;
USES CRT;
VAR
C, Z: ARRAY[1..4] OF REAL;
I: INTEGER;
T, TMAX: REAL;
BEGIN
CLRSCR;
WRITE (‘TMAX=’); READLN (TMAX);
T:=0;
WRITE (‘ВВЕДИТЕ С1, С2, С3, С4’);
READLN(C[1]);
READLN(C[2]);
READLN(C[3]);
READLN(C[4]);
WRITELN (‘T C1 C2 C3 C4’);
WHILE T<=TMAX DO
BEGIN
Z[1]:={ПРАВАЯ ЧАСТЬ ПЕРВОГО УРАВНЕНИЯ СИСТЕМЫ};
Z[2]:={ПРАВАЯ ЧАСТЬ ВТОРОГО УРАВНЕНИЯ СИСТЕМЫ};
Z[3]:={ПРАВАЯ ЧАСТЬ ТРЕТЬЕГО УРАВНЕНИЯ СИСТЕМЫ};
Z[4]:={ПРАВАЯ ЧАСТЬ ЧЕТВЁРТОГО УРАВНЕНИЯ СИСТЕМЫ};
T:=T+0.5;
WRITE (T:5:1);
FOR I:=1 TO 4 DO
WRITE (Z[1]:10:5);
WRITELN;
FOR I:=1 TO 4 DO
C[I]:=Z[I];
END;
READKEY;
END.
Program Eiler;
Uses Crt;
Type Vec=Array[1..6] of Real;
Var y, k, f:Vec;
i, KolKom, KolRec, KolForPrint,count:Byte;
x, xk, h, s:Real;
Fout:text;
NameOut:string[12];
Procedure f_dir(y:vec; var f:vec);
var r:vec;
begin
r[1]:=k[1]*y[1];
r[2]:=k[2]*y[2];
r[3]:=k[3]*y[2];
r[4]:=k[4]*y[3];
r[5]:=k[5]*y[2];
f[1]:=-r[1]+r[5] ;
f[2]:=r[1]-r[2]-r[3]+r[4]-r[5];
f[3]:=r[2]-r[4];
f[4]:=r[3]
end;
BEGIN
ClrScr;
x:=0; xk:=15; h:=0.1;count:=0;s:=0;
write('введите количество компонентов ');readln(KolKom);
For i:=1 to KolKom do
begin
write(введите концентрации компонентов');readln(y[i])
end;
write('введите количество химических реакций ');readln(KolRec);
For i:=1 to KolRec do
begin
write(введите константы химической реакции');readln(k[i])
end;
write('введите число точек для вывода значения на экран');
readln(KolForPrint);
write('введите имя файлов для результатов');
readln(NameOut);
ClrScr;
assign(fout,NameOut);
rewrite(fout);
write('время');write(fout,'время');
For i:=1 to KolKom do begin
write(' y[',i,']');
write(fout,' y[',i,']'); end;
writeln(' сумма '); writeln(fout,' сумма ');
write(x:5:1);write(fout,x:5:1);
for i:=1 to KolKom do
begin s:=s+y[i];write(y[i]:8:5);write(fout,y[i]:8:5) end;
writeln(s:8:2);writeln(fout,s:8:2);
repeat
f_dir(y,f);
For i:=1 to KolKom do y[i]:=y[i]+h*f[i];
x:=x+h;
count:=count+1;
if count=round(xk/h/KolForPrint) then
begin
write(x:5:1);write(fout,x:5:1);s:=0;
for i:=1 to KolKom do
begin
s:=s+y[i];write(y[i]:8:5);write(fout,y[i]:8:5);
end;
writeln(s:8:2);writeln(fout,s:8:2);count:=0
end;
until x>xk;
close(fout);
readkey
end.
ПРОГРАММА МЕТОДА ЭЙЛЕР-КОШИ
Program Eiler_Koshi;
Uses Crt;
Type MAS=Array[1..6] of Real;
Var
y, y0,k, f0,f:MAS;
i, KolKom, KolRec, KolForPrint, count:Byte;
x
xk,
h, s:Real;
Fout:text;
NameOut:string[12];
Procedure f_dir(y:vec; var f:vec);
var r:vec;
begin
r[1]:=k[1]*y[1];
r[2]:=k[2]*y[2];
r[3]:=k[3]*y[2];
r[4]:=k[4]*y[3];
r[5]:=k[5]*y[2];
f[1]:=-r[1]+r[5] ;
f[2]:=r[1]-r[2]-r[3]+r[4]-r[5];
f[3]:=r[2]-r[4];
f[4]:=r[3]
end;
BEGIN
ClrScr;
x:=0; xk:=15; h:=0.1;count:=0;s:=0;
write('введите количество компонентов ');readln(KolKom);
For i:=1 to KolKom do
begin
write(‘введите концентрацию компонента');readln(y[i])
end;
write('введите количество химических реакций ');readln(KolRec);
For i:=1 to KolRec do
begin
write('введите константы химической реакции');readln(k[i])
end;
write('введите число точек для вывода значений на экран ');
readln(KolForPrint);
write('введите имя файла для результатов');
readln(NameOut);
ClrScr;
assign(fout,NameOut);
rewrite(fout);
write(' время');write(fout,' время');
For i:=1 to KolKom do begin
write(' y[',i,']');
write(fout,' y[',i,']'); end;
writeln(' сумма '); writeln(fout,' сумма ');
write(x:5:1);write(fout,x:5:1);
for i:=1 to KolKom do
begin s:=s+y[i];write(y[i]:8:5);write(fout,y[i]:8:5) end;
writeln(s:8:2);writeln(fout,s:8:2);
repeat
f_dir(y,f);
For i:=1 to KolKom do
begin
y0[i]:=y[i];
f0[i]:=f[i];
y[i]:=y0[i]+h*f[i];
end;
f_dir(y,f);
For i:=1 to KolKom do y[i]:=y0[i]+0.5*h*(f0[i]+f[i]);
x:=x+h;
count:=count+1;
if count=round(xk/h/KolForPrint) then
begin
write(x:5:1);write(fout,x:5:1);s:=0;
for i:=1 to KolKom do
begin
s:=s+y[i];write(y[i]:8:5);write(fout,y[i]:8:5);
end;
writeln(s:8:2);writeln(fout,s:8:2);count:=0
end;
until x>xk;
close(fout);
readkey
end.