Смекни!
smekni.com

Решение дифференциального уравнения первого порядка (стр. 2 из 2)

Для этого вторую формулу умножим на

, а третью – на
и сложим с первой. Будем иметь:

Таким образом, с точностью до

имеем приближённую формулу

(3)

Можно показать, что остаточный член формулы (3) равен

где
Аналогично имеем:

и

Отсюда

С другой стороны

Поэтому

Таким образом, с точностью до h5

имеем приближённую формулу

(4)

Можно доказать, что остаточный член формулы (4) есть

где

К формулам (3) и (4) присоединим выражения для производных:

(5)

(6)

Процесс численного дифференцирования уравнения (1) при наличии начального условия (2), использющий формулы (3) и (4), происходит следующим образом. Каким-либо методом вычисляем три начальные строки (начальная таблица):

Из формулы (4) при i=2 получаем первое приближение для

:

(7)

и, пользуясь формулами (5) и (6), находим для соответствующих производных

и
их первые приближения:

и
.

Второе приближение для

определяем при i=2 из формулы (3):

(8)

После этого исправляем значения производных

и
, подсчитывая их вторые приближения:

и
.

Для контроля ещё раз вычисляем по формуле (3) третье приближение

значения
, используя найденные значения
и
.

Если шаг h выбран подходящим, то перещёт не даёт нового результата, и в этом случае можно положить

В противном случае следует уменьшить шаг. Аналогично находятся дальнейшие значения

при i>3.

Для получения начальных значений

и
обычно используют метод последовательных приближений или метод Рунге-Кутта, после чего нужные производные
и
(i=0,1,2) определяются по формулам (5) и (6).

Можна также применить следующий приём: сначала, используя данное начальное значение

, непосредственно вычисляем

и
.

Тем самым будет заполнена первая строка начальной таблицы .

Далее на основании формулы Тейлера приближённо получаем

и, следовательно, можно будит найти

и
.

Пользуясь этими данными, уточняем значение

по формуле (3):

и затем перевычисляем значения

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

Отметим, что если пересчёты элементов строк дают значительные расхождения, то этот приём не является надёжным. В таком случае следует или уменьшить шаг h вычислений, или же обратиться к более точным методам.

В заключение приведём формулы, обеспечивающие более высокую степень точности, но требующие вычисления, кроме второй, ещё и третьей производной искомого решения. А именно, используя Формулу Тейлера и употребляя приём, аналогичный указанному выше, получаем формулы

, (11)

где

, и

, (12)

где

.

Формула (11) употребляется для нахождения первого приближения

; формула (12) даёт уточнённое значение
. Само собой разумеется, что к последним двум формулам целесообразно прибегать тогда, когда форма дифференциального уравнения позволяет сравнительно просто находить вторую и третью производные от искомой функции y.

Приложение

program proizw_w_p;

uses crt;

const epsilon=0.05;

type mas=array[1..100] of real;

nabl=array [1..3] of real;

var i:integer;

x,y,y1,y2:mas;

nabl1,nabl2,nabl3:nabl;

a,h:real;

n:integer;

function f(x, y:real):real;

begin

f:=sqr(x)-sqr(y);

end;

procedure metod(xi, yi, step: real; var rez:real);

var k1, k2, k3, k4:real;

begin

k1:=F(xi,yi);

k2:=F(xi+step/2,yi+k1*step/2);

k3:=F(xi+step/2,yi+k2*step/2);

k4:=F(xi+step,yi+k3*step);

rez:=yi+(step/6)*(k1+2*k2+2*k3+k4)

end;

procedure osn_metod(xi, yi, step:real;var yh22:real;var h:real);

var yh,yh2:real;

begin

repeat

metod(xi, yi,step, yh);

metod(xi, yi, step/2, yh2);

metod(xi, yh2, step/2, yh22);

if abs(yh-yh22)/15>epsilon then step:=h/2;

until abs(yh-yh22)/15<epsilon;

end;

procedure iteraziya(j:integer;xi,h:real);

begin

{первое приближение}

nabl1[1]:=y[j-3]+3*(y[j-1]-y[j-2])+sqr(h)*y2[j-1]-y2[j-2];

{производная первого приближения}

nabl1[2]:=sqr(xi)-sqr(nabl1[1]);

{вторая производная первого приближение}

nabl1[3]:=2*(xi-nabl1[1]*nabl1[2]);

{второе приближение}

nabl2[1]:=y[j-1]+(h/2)*(y1[j-1]+nabl1[2])+((sqr(h))/12)*(nabl1[3]-y2[j-1]);

{производная второго приближения}

nabl2[2]:=sqr(xi)-sqr(nabl2[1]);

{вторая производная второго приближения}

nabl2[3]:=2*(xi-nabl2[1]*nabl2[2]);

{третье приближение}

nabl3[1]:=y[j-1]+(h/2)*(y1[j-1]+nabl2[2])-(sqr(h)/12)*(nabl2[3]-y2[j-1]);

{производная третьего приближения}

nabl3[2]:=sqr(xi)-sqr(nabl3[1]);

{вторая производная третьего приближения}

nabl3[3]:=2*(xi-nabl2[1]*nabl2[2]);

end;

procedure solution(h:real);

begin

{==============Метод Рунге-Кута =================================}

a:=0;

i:=1;

y[1]:=1;

while i<4 do

begin

x[i+1]:=a+i*h;

osn_metod(x[i], y[i], h,y[i+1], h);

inc(i);

end;

{======Окончание метода Рунге-Кута =================================}

{============найдем первые и вторые производные===============}

for i:=1 to 3 do

begin

y1[i]:=sqr(x[i])-sqr(y[i]);

y2[i]:=2*(x[i]-y[i]*y1[i]);

end;

{=================================================================}

for i:=4 to n do

begin

iteraziya(i,x[i],h);

if abs(nabl3[1]-nabl2[1])<epsilon

then

begin

y[i]:=nabl3[1];

y1[i]:=nabl3[2];

y2[i]:=nabl3[3];

end

else

begin

h:=h/2;

if keypressed then halt;

solution(h);

end;

end;

end;

BEGIN

{=====================init==========================================}

clrscr;

write('введите количество значений, которые необходимо вычислить n= ');

readln(n);

h:=0.1;

{==================end of init=========================================}

for i:=1 to n do

begin

x[i]:=(i-1)*h;

end;

solution(h);

for i:=1 to n do

begin

write('y[',i,']= ',y[i],' y"[',i,']= ',y1[i],' y""[',i,']= ',y2[i]);

writeln;

end;

writeln('');

writeln('');

write('Press <enter> to exit....');

readln;

END.