МИНИСТЕРСТВООБРАЗОВАНИЯ И НАУКИ УКРАИНЫ
СУМСКИЙГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
КАФЕДРАИНФОРМАТИКИ
К У Р С О В АЯ Р А Б О Т А
ПО ЧИСЛЕННЫММЕТОДАМ
на тему:
РЕШЕНИЕДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ
ПЕРВОГОПОРЯДКА
Сумы, 2005 г.
1. Метод Адамса
Этот метод численногоинтегрирования разработан Адамсом в 1855г. В последствии этот метод был забыт ивновь открыт в начале века. Популяризация метода Адамса и дальнейшее егоусовершенствование связаны с именем А.Н. Крылова.
Изложим метод Адамсаприменительно к уравнению первого порядка
/> (1)
с начальным условием
/> (2).
Пусть x/>(i=0,1,2,….) – системаравностоящих значений с шагом h и />=/>. Очевидно, имеем
/> (3).
В силу второйинтерполяционной формулы Ньютона с точностью до разностей четвертого порядкаполучаем
/> (4)
где />.
Подставляя выражение (4)в формулу (3) и учитывая, что dx=hdq, будем иметь
/>
Отсюда получаемэкстраполяционную формулу Адамса
/>. (5)
Для начала процесса нужнычетыре начальных значения />, такназываемый начальный отрезок, который определяют исходя из начального условия(2), каким-нибудь численным методом. Можно, например, использовать методРунге-Кутта. Зная эти значения, из уравнения (1) можно найти значенияпроизводных/> и составить таблицуразностей.
/>/>/>/>/>/> (6)
Дальнейшие значения /> (i=4,5,…) искомого решенияможно шаг за шагом вычислять по формуле Адамса, пополняя по мере необходимости таблицу разностей (6).
Для контроля рекомендуетсявычислив первое приближение для /> поформуле
/>
определить />, подсчитать конечныеразности.
/>, />,/> (7)
и затем найти второеприближение по более точной формуле
/> (8)
Если /> и /> отличаются лишь на несколькоединиц последнего сохраняемогодесятичного разряда, то можно положить /> ,а затем, найдя />, перевычисливконечные разности (7). После этого, строго говоря, следует снова найти/> по формуле(8). Поэтому шагh должен быть таким, чтобы этот пересчёт был излишним.
На практике шаг hвыбирают столь малым, чтобы можно было пренебречь членом /> в формуле (8).
Если же расхождениевеличин />и /> значительно, то следуетуменьшить шаг h.
Обычно шаг h уменьшают вдва раза. Покажем, как в этом случае, имея до некоторого значения i таблицувеличин />и/>, />(j/>i) c шагом />, можно просто построитьтаблицу величин />(k=0,1,2…) сшагом />. Для кратности введениясокращенные обозначения:
/>(k=0,1,2…).
На основе формулы (4)будем иметь
/>, (9)
где />. Отсюда, полагая j=i-2 иq=1/2 и учитывая, что />,находим
/>. (10)
Аналогично при j=i-1,q=1/2 из формулы (9) получаем, что аргументу /> соответствует значение
/>. (11)
Что касается значений /> и /> , то они имеются в старойтаблице. После этого составляем начальный отрезок для новой таблицы./> и находим конечныеразности:
/> (k=-3,-2,-1),
/> (k=-3,-2),
/> (k=-3,).
Дальше таблицапродолжается обычным путём, посредством соответствующей модификации формулы(5):
/>,
/> (j=0,1,2,…).
Для работы на электронныхсчётчиках машинах формулу Адамса (5) выгодно применять в раскрытом виде.Учитывая, что
/>
/>
/>
после приведения подобныхчленов имеем
/>,
причём />.
2. Методы, основанные на применениипроизводных высших порядков
До сих пор для численногоинтегрирования дифференциального уравнения первого порядка
/> (1)
с начальным условием
/>(2)
мы применяли формулы, вкоторых явно используется лишь первая производная /> искомогорешения.
Однако если использоватьформулы, явно содержащие производные высших порядков от искомого решения, томожно указать методы, дающие более точный результат на данном промежутке безувеличения числа шагов.
Выведем соответствующиеформулы, предполагая, что правая часть уравнения (1) дифференцируема достаточное число раз.
Пусть /> — значения искомого решенияy=y(x) и, соответственно, значения его производных первого и второго порядков вточках />. Располагая величины
/>
в ряды по степеням h,находим:
/>
/>
/>
Из полученных формулисключим члены, содержащие />/> и />.
Для этого вторую формулуумножим на />, а третью – на /> и сложим с первой. Будемиметь:
/>
Таким образом, сточностью до /> имеем приближённую формулу
/> (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] ofreal;
nabl=array [1..3] ofreal;
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);
ifabs(yh-yh22)/15>epsilon then step:=h/2;
untilabs(yh-yh22)/15
end;
procedureiteraziya(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;
proceduresolution(h:real);
begin
{==============МетодРунге-Кута =================================}
a:=0;
i:=1;
y[1]:=1;
while i
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);
ifabs(nabl3[1]-nabl2[1])
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 ofinit=========================================}
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 to exit....');
readln;
END.