Численный расчет дифференциальных уравнений

Міністерство освіти України

ДАЛПУ

Кафедра автоматизації

технологічних процесів і приладобудування

КУРСОВА РОБОТА

з курсу “Математичне моделювання на ЕОМ”

на тему “Розв’язок диференціального рівняння

виду а_>п>у^(п)+а_>п-1>у^(п-1)+…+а_>1>у¹+а_>0>у=кх при заданих

початкових умовах з автоматичним вибором кроку

методом Ейлера”

Виконала студентка групи БА-4-97

Богданова Ольга Олександрівна

Холоденко Вероніка Миколаївна

Перевірила Заргун Валентина Василівна

1998

Блок-схема алгоритма

Блок-схема алгоритма

начало

у^/=f(x,y)

y(x_>0>)=y_>0 >

x_>0>, x_>0>+a

h, h/2

k:=0

x_>k+1/2>:=x_>k>+h/2

y_>k+1/2>:=y_>k>+f(x_>k, >y_>k>)h/2

_> >α_>k>:=_> >f(x_>k+1/2,> y_>k+1/2>)

x_>k+1>:=x_>k>+h

y_>k+1>:=y_>k>+α_>k>h

_>нет> k:=n

_{> да}>

_> >

x_>0>, y_>0>,

x1, y1…

x_>n>,_> >y_>n>

_> > конец

_> >

ПОСТАНОВКА ЗАДАЧИ И МЕТОД РЕШЕНИЯ

Решить дифференциальное уравнение у^/=f(x,y) численным методом - это значит для заданной последовательности аргументов х_>0>, х_>1>…, х_>n> и числа у_>0>, не определяя функцию у=F(x), найти такие значения у_>1>, у_>2>,…, у_>n>, что у_>i>=F(x_>i>)(i=1,2,…, n) и F(x_>0>)=y_>0>.

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

У=F(x) получить таблицу значений этой функции для заданной последовательности аргументов. Величина h=x_>k>-x_>k>_>-1> называется шагом интегрирования.

Метод Эйлера относиться к численным методам, дающим решение в виде таблицы приближенных значений искомой функции у(х). Он является сравнительно грубым и применяется в основном для ориентировочных расчетов. Однако идеи, положенные в основу метода Эйлера, являются исходными для ряда других методов.

Рассмотрим дифференциальное уравнение первого порядка

y^/=f(x,y) (1)

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

x=x_>0>, y(x_>0>)=y_>0 > (2)

Требуется найти решение уравнения (1) на отрезке [а,b].

Разобьем отрезок [a, b] на n равных частей и получим последовательность х_>0>, х_>1>, х_>2>,…, х_>n>, где x_>i>=x_>0>+ih (i=0,1,…, n), а h=(b-a)/n-шаг интегрирования.

В методе Эйлера приближенные значения у(х_>i>)»y_>i> вычисляются последовательно по формулам у_>i>+hf(x_>i>, y_>i>) (i=0,1,2…).

При этом искомая интегральная кривая у=у(х), проходящая через точку М_>0>(х_>0>, у_>0>), заменяется ломаной М_>0>М_>1>М_>2>… с вершинами М_>i>(x_>i>, y_>i>) (i=0,1,2,…); каждое звено М_>i>M_>i+1> этой ломаной, называемой ломаной Эйлера, имеет направление, совпадающее с направлением той интегральной кривой уравнения (1), которая проходит через точку М_>i>.

Если правая часть уравнения (1) в некотором прямоугольнике R{|x-x_>0>|£a, |y-y_>0>|£b}удовлетворяет условиям:

|f(x, y_>1>)- f(x, y_>2>)| £ N|y_>1>-y_>2>| (N=const),

|df/dx|=|df/dx+f(df/dy)| £ M (M=const),

то имеет место следующая оценка погрешности:

|y(x_>n>)-y_>n>| £ hM/2N[(1+hN)ⁿ-1], (3)

где у(х_>n>)-значение точного решения уравнения(1) при х=х_>n>, а у_>n>- приближенное значение, полученное на n-ом шаге.

Формула (3) имеет в основном теоретическое применение. На практике иногда оказывается более удобным двойной просчет: сначала расчет ведется с шагом h, затем шаг дробят и повторный расчет ведется с шагом h/2. Погрешность более точного значения у_>n>^* оценивается формулой

|y_>n>-y(x_>n>)|»|y_>n>^*-y_>n>|.

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

Модифицированный метод Эйлера более точен.

Рассмотрим дифференциальное уравнение (1) y^/=f(x,y)

с начальным условием y(x_>0>)=y_>0>. Разобьем наш участок интегрирования на n

равных частей. На малом участке [x_>0>,x_>0>+h]

_>у >интегральную кривую заменим прямой

_> >N_>k>^/ _>y=y(x) > линией. Получаем точку М_>к>(х_>к>,у_>к>).

_> >

М_>к> М_>к>^/

y_>k+1 >

y_>k>

х_>к> х_>к1>_>/2> x_>k+h=>x_>k1> х

Через М_>к >проводим касательную: у=у_>к>=f(x_>k>,y_>k>)(x-x_>k>).

Делим отрезок (х_>к>,х_>к1>) пополам:

x_>Nk>^/=x_>k>+h/2=x_>k+1/2>

_> > y_>Nk>^/=y_>k>+f(x_>k>,y_>k>)h/2=y_>k>+y_>k+1/2>

Получаем точку N_>k>^/. В этой точке строим следующую касательную:

y(x_>k+1/2>)=f(x_>k+1/2>, y_>k+1/2>)=α_>k>

Из точки М_>к> проводим прямую с угловым коэффициентом α_>к> и определяем точку пересечения этой прямой с прямой Х_>к1>. Получаем точку М_>к>^/. В качестве у_>к+1> принимаем ординату точки М_>к>^/. Тогда:

у_>к+1>=у_>к>+α_>к>h

x_>k+1>=x_>k>+h

(4) α_>k>=f(x_>k+h/2>, y_>k>+f(x_>k>,Y_>k>)h/2)

y_>k>=y_>k-1>+f(x_>k-1>,y_>k-1>)h

(4)-рекурентные формулы метода Эйлера.

Сначала вычисляют вспомогательные значения искомой функции у_>к+1>_>/>_>2> в точках х_>к+1>_>/2>, затем находят значение правой части уравнения (1) в средней точке y^/_>k>_>+1>_>/2>=f(x_>k+1/2>, y_>k+1/2>) и определяют у_>к+1>.

Для оценки погрешности в точке х_>к> проводят вычисления у_>к> с шагом h, затем с шагом 2h и берут 1/3 разницы этих значений:

| у_>к>^*-у(х_>к>)|=1/3(y_>k>^*-y_>k>),

где у(х)-точное решение дифференциального уравнения.

Таким образом, методом Эйлера можно решать уравнения любых порядков. Например, чтобы решить уравнение второго порядка y^//=f(y^/,y,x) c начальными условиями y^/(x_>0>)=y^/_>0>, y(x_>0>)=y_>0>, выполняется замена:

y^/=z

z^/=f(x,y,z)

Тем самым преобразуются начальные условия: y(x_>0>)=y_>0>, z(x_>0>)=z_>0>, z_>0>=y^/_>0>.

РЕШЕНИЕ КОНТРОЛЬНОГО ПРИМЕРА

Приведем расчет дифференциального уравнения первого, второго и третьего порядка методом Эйлера

1. Пусть дано дифференциальное уравнение первого порядка:

y^/=2x-y

Требуется найти решение на отрезке [0,1] c шагом h=(1-0)/5=0,2

Начальные условия: у_>0>=1;

Пользуясь рекурентными формулами (4), находим:

1). x_>1>=0,2; х_>1>_>/2>=0,1; y(x_>1>)=y(x_>0>)+α_>0>h; y(x_>1/2>)=y(x_>0>)+f(x_>0>,y_>0>)h/2;

f(x_>0>,y_>0>)=2*0-1=-1

y(x_>1/2>)=1-1*0,1=0,9

α_>0>=2*0,1-0,9=-0,7

y_>1>=1-0,1*0,2=0,86

2). y(x_>2>)=y(x_>1>)+α_>1>h; x_>2>=0,2+0,2=0,4; x_>1+1/2>=x_>1>+h/2=0,2+0,1=0,3

y(x_>1+1/2>)=y(x_>1>)+f(x_>1>,y(x_>1>))h/2

f(x_>1>,y_>1>)=2*0,2-0,86=-0,46

y(x_>1+1/2>)=0,86-0,46*0,1=0,814

α_>1>=2*0,3-0,814=-0,214

y_>2>=0,86-0,214*0,2=0,8172

3). x_>3>=0,4+0,2=0,6; x_>2+1/2>=x_>2>+h/2=0,4+0,1=0,5

f(x_>2>,y_>2>)=2*0,4-0,8172=-0,0172

y_>2+1/2>=0,8172-0,0172*0,1=0,81548

α_>2>=2*0,5-0,81548=0,18452

y_>3>=0,8172+0,18452*0,2=0,854104

4).x_>4>=0,8; x_>3+1/2>=x_>3>+h/2=0,6+0,1=0,7

f(x_>3>,y_>3>)=2*0,6-0,854104=0,345896

y_>3+1/2>=0,854104+0,345896*0,1=0,8886936

α_>3>=2*0,7-0,89=0,5113064

y_>4>=0,854104+0,5113064*0,2=0,95636528

5).x_>5>=1; x_>4+1/2>=0,8+0,1=0,9

f(x_>4>,y_>4>)=2*0,8-0,956=0,64363472

y_>4+1/2>=0,956+0,643*0,1=1,020728752;

α_>4>=2*0,9-1,02=0,779271248

y_>5>=0,956+0,7792*0,2=1,11221953

2. Дано уравнение второго порядка:

y^//=2x-y+y^/

Находим решение на том же отрезке [0,1] c шагом h=0,2;

Замена: y^/=z

z^/=2x-y+z

Начальные условия: у_>0>=1

z_>0>=1

1).x_>1>=0,2; x_>1/2>=0,1

y(z_>1>)=y(z_>0>)+α_>0>h z(x_>1>,y_>1>)=z(x_>0>,y_>0>)+β_>0>h

y(z_>1/2>)=y(z_>0>)+f(z_>0>,y_>0>)h/2 z(x_>1/2>,y_>1/2>)=z(x_>0>,y_>0>)+f(x_>0>,y_>0>,z_>0>)h/2

f(z_>0>,y_>0>)=f_>10>=1 f(x_>0>,y_>0>,z_>0>)=f_>20>=2*0-1+1=0

y_>1/2>=1+1*0,1=1,1 z_>1/2>=1+0*0,1=1

α_>0>=z_>0>=1 β_>0>=2*0,1-1,1+1=0,1

y_>1>=1+0,2*1=1,2 z_>1>=1+0,2*0,1=1,02

2).x_>2>+0,4; x_>1+1/2>=0,3

f_>11>=z_>1>=1,02 f_>21>=2*0,2-1,2+1,02=0,22

y_>1+1/2>=1,2+1,02*0,1=1,1 z_>1+1/2>=1,02+0,22*0,1=1,042

α_>1>=z_>1+1/2>=1,042 β_>1>=2*0,3-1,302+1,042=0,34

y_>2>=1,2+1,042*0,2=1,4084 z_>2>=1.02+0,34*0,2=1,088

3).x_>3>=0,6; x_>2+1/2>=0,5

f_>12>=z_>2>=1,088 f_>22>=2*0,4-1,4084+1,088=0,4796

y_>2+1/2>=1,4084+1,088*0,1=1,5172 z_>2+1/2>=1,088+0,4796*0,1=1,13596

α_>2>=z_>2+1/2>=1,13596 β_>2>=2*0,5-1,5172+1,13596=0,61876

y_>3>=1,4084+1,136*0,2=1,635592 z_>3>=1,088+0,61876*0,2=1,211752

4).x_>4>=0,8; x_>3+1/2>=0,7

f_>13>=z_>3>=1,211752 f_>23>=2*0,6-1,636+1,212=0,77616

y_>3+1/2>=1,636+1,212*0,1=1,7567672 z_>3+1/2>=1,212+0,776*0,1=1,289368

α_>3>=z_>3+1/2>=1,289368 β_>3>=2*0,7-1,7568+1,289=0,9326008

y_>4>=1,6+1,289*0,2=1,8934656 z_>4>=1,212+0,93*0,2=1,39827216

5).x_>5>=1; y_>4+1/2>=0,9

f_>14>=z_>4>=1,39827216 f_>24>=2*0,8-1,893+1,398=1,10480656

y_>4+1/2>=1,893+1,398*0,1=2,0332928 z_>4+1/2>=1,398+1,105*0,1=1,508752816

α_>4>=z_>4+1/2>=1,508752816 β_>4>=2*0,9-2,03+1,5=1,27546

y_>5>=1,893+1,5*0,2=2,195216163 z_>5>=1,398+1,275*0,2=1,65336416

3. Чтобы решить уравнение третьего порядка

y^///=2x-y-y^/+y^//

на отрезке [0,1], с шагом h=0,2 и начальными условиями

y_>0>^//=1

y_>0>^/=1

y_>0>=1

необходимо сделать 3 замены: y^/=a y_>0>^/=a_>0>=1

y^//=a^/=b y_>0>^//=b_>0>=1

b^/=2x-y-a+b

1).x_>1>=0,2; x_>1/2>=0,1

y(a_>1>)=y(a_>0>)+a_>0>h y(a_>1/2>)=y(a_>0>)+f_>10>h/2

a(b_>1>)=a(b_>0>)+β_>0>h a(b_>1/2>)=a(b_>0>)+f_>20>h/2

b(x_>1>,y_>1>,a_>1>)=b(x_>0>,y_>0>,a_>0>)+γ_>0>h b(x_>1/2>,y_>1/2>,a_>1/2>)=b(x_>0>,y_>0>,a_>0>)+f_>30>h/2

f_>10>=f(a_>0>,y(a_>0>))=1 y_>1/2>=1+1*0,1=1,1

f_>20>=f(b_>0>,a(b_>0>))=1 a_>1/2>=1+1*0,1=1,1

f_>30>=f(x_>0>,y_>0>,a_>0>,b_>0>)=-1 b_>1/2>=1-1*0,1=0,9

α_>0>=a_>1/2>=1,1 y(a_>1>)=1+1,1*0,2=1,22

β_>0>=b_>1/2>=0,9 a(b_>1>)=1+0,9*0,2=1,18

γ_>0>=2*0,1-1,1-1,1+0,9=-1,1 b(x_>1>,y_>1>,a_>1>)=1-1,1*0,2=0,78

2).x_>2>=0,4; x_>1+1/2>=x_>1>+h/2=0,3

f_>11>=a_>1>=1,18 y_>1+1/2>=1,22+1,18*0,1=1.338

f_>21>=b_>1>=0,78 a_>1+1/2>=1,18+0,78*0,1=1,258

f_>31>=2*0,2-1,22-1,18+0,78=-1,22 b_>1+1/2>=-1,22*0,1+0,78=0,658

α_>1>=a_>1+1/2>=1,258 y_>2>=1,22+1,258*0,2=1,4716

β_>1>=b_>1+1/2>=0,658 a_>2>=1,18+0,658*0,2=1,3116

γ_>1>=2*0,3-1,338-1,258+0,658=-1,338 b_>2>=0,78-1,338*0,2=0,5124

3).x_>3>=0,6; x_>2+1/2>=0,5

f_>12>=a_>2>=1,3116 y_>2+1/2>=1,47+1,3*0,1=1,60276

f_>22>=b_>2>=0,5124 a_>2+1/2>=1,3116+0,5*0,1=1.36284

f_>32>=2*0,4-1,47-1,31+0,512=-1,4708 b_>2+1/2>=0,4-1,4*0,1=0,36542

α_>2>=1,36284 y_>3>=1,4716+1,3116*0,2=1,744168

β_>2>=0,36542 a_>3>=1,3116+0,3654*0,2=1,384664

γ_>2>=2*0,5-1,6-1,36+0,365=-1,60018 b_>3>= 0,51-1,60018*0,2=0,192364

4).x_>4>=0,8; x_>3+1/2>=0,7

f_>13>=1,384664 y_>3+1/2>=1,74+1,38*0,1=1,8826364

f_>23>=0,192364 a_>3+1/2>=1,38+0,19*0,1=1,4039204

f_>33>=2*0,6-1,7-1,38+0,19=-1,736488 b_>3+1/2>=0,19-1,7*0,1=0,0187152

α_>3>=1,4039204 y_>4>=1,74+1,4*0,2=2,0249477

β_>3>=0,0187152 a_>4>=1,38+0,9187*0,2=1,388403

γ_>3>=2*0,7-1,88-1,4+0,0187=-1,8678416 b_>4>=0,192-1,87*0,2=-0,1812235

5).x_>4>=1; x_>4+1/2>=0,9

f_>14>=1,388403 y_>4+1/2>=2,02+1,388*0,1=2,16379478

_> > f_>24>=-0,1812235 a_>4+1/2>=1,4-0.181*0,1=1,370306608

f_>34>=2*0,8-2,02-1,388-0,18=-1,9945834 b_>4+1/2>=-0,18-1,99*0,1=-0,38066266

α_>4>=1,3703 y_>5>=2,02+1,37*0,2=2,2990038

β_>4>=-0,38066 a_>5>=1,388-0,38*0,2=1,3122669

γ_>4>=2*0,9-2,16-1,37-0,38=-2,114764056 b_>5>=-0,181-2,1*0,2=-0,6041734

Программа на Turbo Pascal

uses crt,pram,kurs1_1;

var

yx,xy,l,v,p,ff,ay,by,x:array [0..10] of real;

y,a,b:array[0..10,0..1] of real;

i,n,o:integer;

c,d,h,k:real;

label

lap1;

begin

screen1;

clrscr;

writeln('введите наивысший порядок производной не больше трех ');

readln(n);

if n=0 then begin

writeln('это прямолинейная зависимость и решается без метода Эйлера ');

goto lap1;end;

writeln('введите коэффициенты {a0,a1}');

for i:=0 to n do

readln(l[i]);

if (n=1) and (l[1]=0) or (n=2) and (l[2]=0) or (n=3) and (l[3]=0) then begin

writeln('деление на ноль');

goto lap1;

end;

writeln('введите коэффициент при x');

readln(k);

writeln('введите отрезок ');

readln(c,d);

o:=5;

h:=abs(d-c)/o;

writeln('шаг=',h:1:1);

writeln('задайте начальные условия y(x)= ');

for i:=0 to n-1 do

readln(v[i]);

if n=3 then begin

yx[0]:=v[0];

ay[0]:=v[1];

by[0]:=v[2];

p[0]:=(k*c-l[0]*v[0]-l[1]*v[1]-l[2]*v[2])/l[3];

x[0]:=c;

gotoxy(32,1);

write(' ');

gotoxy(32,2);

write(' x y a b ');

gotoxy(32,3);

write(' ',c:7:7,' ',yx[0]:7:7,' ',ay[0]:7:7,' ',by[0]:7:7,' ');

for i:=0 to o-1 do begin

x[i]:=x[i]+h/2;

y[i,1]:=yx[i]+(h/2)*ay[i];

a[i,1]:=ay[i]+(h/2)*by[i];

b[i,1]:=by[i]+(h/2)*p[i];

ff[i]:=(k*x[i]-l[0]*y[i,1]-l[1]*a[i,1]-l[2]*b[i,1])/l[3];

xy[i]:=x[i]+h/2;

yx[i+1]:=yx[i]+h*a[i,1];

ay[i+1]:=ay[i]+h*b[i,1];

by[i+1]:=by[i]+h*ff[i];

x[i+1]:=x[i]+h/2;

p[i+1]:=(k*xy[i]-l[0]*yx[i+1]-l[1]*ay[i+1]-l[2]*by[i+1])/l[3];

end;

for i:=0 to o-1 do begin

gotoxy(32,4+i);

write(' ',xy[i]:7:7,' ',yx[i+1]:7:7,' ',ay[i+1]:7:7,' ',by[i+1]:7:7,' ');

end;

gotoxy(32,4+o);

write(' ');

end;

if n=2 then begin

x[0]:=c;

yx[0]:=v[0];

ay[0]:=v[1];

p[0]:=(k*c-l[0]*yx[0]-l[1]*v[1])/l[2];

gotoxy(32,1);

write(' ');

gotoxy(32,2);

write(' x y a ');

gotoxy(32,3);

write(' ',c:7:7,' ',yx[0]:7:7,' ',ay[0]:7:7,' ');

for i:=0 to o-1 do begin

x[i]:=x[i]+h/2;

y[i,1]:=yx[i]+(h/2)*ay[i];

a[i,1]:=ay[i]+(h/2)*p[i];

ff[i]:=(k*x[i]-l[0]*y[i,1]-l[1]*a[i,1])/l[2];

xy[i]:=x[i]+h/2;

yx[i+1]:=yx[i]+h*a[i,1];

ay[i+1]:=ay[i]+h*ff[i];

x[i+1]:=x[i]+h/2;

p[i+1]:=(k*xy[i]-l[0]*yx[i+1]-l[1]*ay[i+1])/l[2];

end;

for i:=0 to o-1 do begin

gotoxy(32,4+i);

write(' ',xy[i]:7:7,' ',yx[i+1]:7:7,' ',ay[I+1]:7:7,' ');

end;

gotoxy(32,4+o);

write(' ');

end;

if n=1 then begin

x[0]:=c;

yx[0]:=v[0];

p[0]:=(k*x[0]-l[0]*yx[0])/l[1];

for i:=0 to o-1 do begin

x[i]:=x[i]+h/2;

y[i,1]:=yx[i]+(h/2)*p[i];

xy[i]:=x[i]+h/2;

ff[i]:=(k*x[i]-l[0]*y[i,1])/l[1];

yx[i+1]:=yx[i]+h*ff[i];

x[i+1]:=x[i]+h/2;

p[i+1]:=(k*xy[i]-l[0]*yx[i+1])/l[1];

end;

gotoxy(32,1);

write(' ');

gotoxy(32,2);

write(' x y ');

gotoxy(32,3);

write(' ',c:7:7,' ',yx[0]:7:7,' ');

for i:=0 to o-1 do begin

gotoxy(32,4+i);

write(' ',xy[i]:7:7,' ',yx[i+1]:7:7,' ');

end;

gotoxy(32,o+4);

write(' ');

end;

lap1:readln;

pramo;

delay(10000);

clrscr;

end.



ЗАПУСК ПРОГРАММЫ НА ВЫПОЛНЕНИЕ

Программа находится в файле kursova1.pas, и имеет 2 модуля, в которых содержатся заставки. Модули находятся в файлах pram.tpu и kurs1_1.tpu.

Для запуска файла kursova1.pas в Turbo Pascal необходимо нажать F9. Появится первая заставка, далее нажать enter и ввести все необходимые начальные условия: порядок производной, коэффициенты при членах рада, отрезок и начальные значения у(х_>0>). На экране выводится шаг вычисления и таблица с ответами. После нажатия enter выводится вторая заставка, после чего мы возвращаемся к тексту программы.

ОПИСАНИЕ ПРОГРАММЫ

1 – ввод данных, используемых в программе

2 – использование метки, очистка экрана, ввод требований, решение

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

условий

3 – присвоение начальных условий для дифференциального уравнения

третьего порядка

4 – вывод таблицы со значениями

5 – ввод формул метода Эйлера для уравнения третьего порядка

6 – присвоение начальных условий для решения дифференциального

уравнения второго порядка

7 – вывод таблицы для уравнения второго порядка

8 – формулы метода Эйлера для уравнения второго порядка

9 – начальные условия для дифференциального уравнения первого порядка

10 – формулы метода Эйлера для решения уравнения первого порядка

11 – вывод таблицы

12 – обращение к метке, задержка для просмотра результатов, очистка

экрана, конец программы.