3.1 Метод эйлера.
Пусть дано дифференциальное уравнение y'=f(x,y) (1), с начальными условиями y(x0)=y0 .
Пусть y=y(x) – искомое точное решение. Интегральная кривая проходит через точку (x0,y0).
Найдем приближенные значения функции в точках x1, x2, …, xn. Построим систему равноотстоящих точек x0, x1=x0+h, …, xn=b. Проведём прямые x=x0, x=x1, …, x=b.
y2
y1
y0
Рассмотрим отрезок [x0,x1]. На этом отрезке есть одна точка, которая принадлежит искомой кривой – это точка А (x0,y0). Заменим дугу искомой кривой y=y(x) на отрезке [x0,x1] касательной к ней, проведенной в точке (x0,y0). В качестве y1 возьмём ординату точки пересечения прямой x=x1 к касательной.
Очевидно y1=y0+h•tgα0. Но tgα0=y'(x0), т.е. y1=y0+h•y'(x0). Но из уравнения (1) следует, что
y'(x0)=f(x0,y0).
Итак, получаем y1=y0+h•f(x0,y0), x1=x0+h .
Предположим теперь, что точка (x1,y1) принадлежит искомой кривой. В этой точке опять проведём касательную к графику функции до пересечения с прямой x=x2 .
Тогда аналогично:
y2= y1+ h f(x1,y1);
x2= x1+ h .
Продолжая и так далее, получим систему значений y1,y2,…,yn, которые и будут приближенными значениями функции y=y(x) в точках x1,x2,…xn.
Итак, расчётные формулы метода Эйлера:
ym+1= ym+ h f(xm,ym);
xm+1= xm+ h
Для системы дифференциальных уравнений
y'I = fi(x,y1,y2,…yk)
yi(x0) = yi0 i = 1,…,K
расчетные формулы записываются аналогично
yim+1 = yim + hfi(x,y1m,y2m,…ykm);
xm+1 = xm+ h,
здесь i – номер уравнения в системе, n – номер шага.
Метод Эйлера является грубым методом, ошибка, которую мы допускаем, на каждом шаге, пропорциональна h, т.е. |y(xk) - yk|≈ h.
Чтобы повысить точность вычислений, используют некоторые усовершенствованные методы.
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.