Математика

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



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


                                 ДАЛПУ



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

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



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



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

        на тему “Розв’язок диференціального рівняння
             виду  апу(п)+ап-1у(п-1)+…+а1у1+а0у=кх  при заданих
               початкових умовах з автоматичним вибором кроку
                               методом Ейлера”



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


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

                                                          Холоденко
Вероніка Миколаївна
                                                          Перевірила Заргун
Валентина Василівна



                                    1998


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



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

                                   начало

                                  у/=f(x,y)
                                  y(x0)=y0
                                              x0, x0+a


                                   h, h/2


                                    k:=0


                                        xk+1/2:=xk+h/2
                               yk+1/2:=yk+f(xk, yk)h/2
                                                             ?k:= f(xk+1/2,
yk+1/2)
                                        xk+1:=xk+h
                                         yk+1:=yk+?kh


                                    нет         k:=n

        да

                                                x0, y0,
                                   x1, y1…
                                                xn, yn


                                                конец



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


     Решить дифференциальное уравнение у/=f(x,y)  численным  методом  -  это
значит для заданной последовательности аргументов х0, х1…, хn  и  числа  у0,
не определяя  функцию  у=F(x),  найти  такие  значения  у1,  у2,…,  уn,  что
уi=F(xi)(i=1,2,…, n) и F(x0)=y0.
Таким образом, численные методы позволяют вместо нахождения функции
У=F(x)   получить   таблицу   значений    этой    функции    для    заданной
последовательности   аргументов.   Величина   h=xk-xk-1   называется   шагом
интегрирования.
      Метод Эйлера относиться к численным методам,  дающим  решение  в  виде
таблицы  приближенных  значений   искомой   функции   у(х).   Он    является
сравнительно грубым и применяется в основном для  ориентировочных  расчетов.
Однако идеи, положенные в основу метода Эйлера, являются исходными для  ряда
других методов.

      Рассмотрим дифференциальное уравнение первого порядка
                                       y/=f(x,y)                (1)
с начальным условием
                                     x=x0, y(x0)=y0              (2)
Требуется найти решение уравнения (1) на отрезке [а,b].
Разобьем отрезок [a, b] на n равных частей и получим последовательность  х0,
х1, х2,…, хn, где xi=x0+ih (i=0,1,…, n), а h=(b-a)/n-шаг интегрирования.
        В  методе  Эйлера   приближенные   значения   у(хi)(yi   вычисляются
последовательно по формулам уi+hf(xi, yi) (i=0,1,2…).
При этом искомая интегральная кривая у=у(х), проходящая через  точку  М0(х0,
у0), заменяется ломаной М0М1М2… с вершинами Мi(xi, yi)  (i=0,1,2,…);  каждое
звено МiMi+1 этой ломаной, называемой  ломаной  Эйлера,  имеет  направление,
совпадающее с направлением той интегральной кривой  уравнения  (1),  которая
проходит через точку Мi.
     Если правая часть уравнения (1) в некотором прямоугольнике R{|x-x0|(a,
|y-y0|(b}удовлетворяет условиям:

                           |f(x, y1)- f(x, y2)| ( N|y1-y2|  (N=const),
                           |df/dx|=|df/dx+f(df/dy)| ( M  (M=const),

то имеет место следующая оценка погрешности:
                          |y(xn)-yn| ( hM/2N[(1+hN)n-1],            (3)

где  у(хn)-значение  точного  решения   уравнения(1)   при   х=хn,   а   уn-
приближенное значение, полученное на n-ом шаге.
Формула (3) имеет в основном теоретическое применение.  На  практике  иногда
оказывается более удобным двойной просчет: сначала расчет  ведется  с  шагом
h, затем шаг дробят и повторный расчет ведется  с  шагом   h/2.  Погрешность
более точного значения уn* оценивается формулой
                                          |yn-y(xn)|(|yn*-yn|.

    Метод  Эйлера  легко  распространяется   на   системы   дифференциальных
уравнений и на дифференциальные уравнения высших порядков. Последние  должны
быть предварительно приведены к системе дифференциальных  уравнений  первого
порядка.
     Модифицированный метод Эйлера более точен.
Рассмотрим дифференциальное уравнение (1)  y/=f(x,y)
с начальным условием y(x0)=y0. Разобьем  наш  участок  интегрирования  на  n

                                          равных частей.  На  малом  участке
[x0,x0+h]
у
интегральную кривую заменим прямой
                     Nk/                            y=y(x)         линией.
Получаем точку Мк(хк,ук).


           Мк          Мк/
                        yk+1
            yk

             хк хк1/2 xk+h=xk1        х


  Через Мк проводим касательную:  у=ук=f(xk,yk)(x-xk).
Делим отрезок (хк,хк1) пополам:
                                          xNk/=xk+h/2=xk+1/2

yNk/=yk+f(xk,yk)h/2=yk+yk+1/2
Получаем точку Nk/. В этой точке строим следующую касательную:
                                           y(xk+1/2)=f(xk+1/2, yk+1/2)=?k
Из точки Мк проводим прямую с угловым коэффициентом ?к  и  определяем  точку
пересечения этой прямой с прямой Хк1. Получаем точку Мк/.  В  качестве  ук+1
принимаем ординату точки Мк/. Тогда:
                                           ук+1=ук+?кh
                                           xk+1=xk+h
                          (4)            ?k=f(xk+h/2, yk+f(xk,Yk)h/2)
                                           yk=yk-1+f(xk-1,yk-1)h
(4)-рекурентные формулы метода Эйлера.
      Сначала вычисляют вспомогательные значения искомой  функции  ук+1/2  в
точках хк+1/2, затем находят значение правой части уравнения (1)  в  средней
точке y/k+1/2=f(xk+1/2, yk+1/2) и определяют ук+1.
      Для оценки погрешности в точке хк проводят вычисления ук с шагом h,
затем с шагом 2h и берут 1/3 разницы этих значений:
                                          | ук*-у(хк)|=1/3(yk*-yk),
где у(х)-точное решение дифференциального уравнения.


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


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


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

1.  Пусть дано дифференциальное уравнение первого порядка:
                                   y/=2x-y
Требуется найти решение на отрезке [0,1] c шагом h=(1-0)/5=0,2
Начальные условия: у0=1;
Пользуясь рекурентными формулами (4), находим:
1). x1=0,2;  х1/2=0,1;     y(x1)=y(x0)+?0h;   y(x1/2)=y(x0)+f(x0,y0)h/2;
      f(x0,y0)=2(0-1=-1
      y(x1/2)=1-1(0,1=0,9
      ?0=2(0,1-0,9=-0,7
       y1=1-0,1(0,2=0,86

2). y(x2)=y(x1)+?1h;   x2=0,2+0,2=0,4;   x1+1/2=x1+h/2=0,2+0,1=0,3
     y(x1+1/2)=y(x1)+f(x1,y(x1))h/2
     f(x1,y1)=2(0,2-0,86=-0,46
     y(x1+1/2)=0,86-0,46(0,1=0,814
     ?1=2*0,3-0,814=-0,214
     y2=0,86-0,214*0,2=0,8172

3). x3=0,4+0,2=0,6;   x2+1/2=x2+h/2=0,4+0,1=0,5
     f(x2,y2)=2*0,4-0,8172=-0,0172
     y2+1/2=0,8172-0,0172*0,1=0,81548
     ?2=2*0,5-0,81548=0,18452
     y3=0,8172+0,18452*0,2=0,854104

 4).x4=0,8;   x3+1/2=x3+h/2=0,6+0,1=0,7
     f(x3,y3)=2*0,6-0,854104=0,345896
     y3+1/2=0,854104+0,345896*0,1=0,8886936
     ?3=2*0,7-0,89=0,5113064
     y4=0,854104+0,5113064*0,2=0,95636528

5).x5=1;   x4+1/2=0,8+0,1=0,9
    f(x4,y4)=2*0,8-0,956=0,64363472
    y4+1/2=0,956+0,643*0,1=1,020728752;
    ?4=2*0,9-1,02=0,779271248
    y5=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
                                           z0=1

1).x1=0,2;   x1/2=0,1
   y(z1)=y(z0)+?0h                                   z(x1,y1)=z(x0,y0)+?0h

   y(z1/2)=y(z0)+f(z0,y0)h/2
z(x1/2,y1/2)=z(x0,y0)+f(x0,y0,z0)h/2
   f(z0,y0)=f10=1                                       f(x0,y0,z0)=f20=2*0-
1+1=0
   y1/2=1+1*0,1=1,1                                 z1/2=1+0*0,1=1
   ?0=z0=1                                                ?0=2*0,1-
1,1+1=0,1
   y1=1+0,2*1=1,2                                   z1=1+0,2*0,1=1,02

2).x2+0,4;   x1+1/2=0,3
    f11=z1=1,02                                         f21=2*0,2-
1,2+1,02=0,22
    y1+1/2=1,2+1,02*0,1=1,1                     z1+1/2=1,02+0,22*0,1=1,042
    ?1=z1+1/2=1,042                                   ?1=2*0,3-
1,302+1,042=0,34
    y2=1,2+1,042*0,2=1,4084                  z2=1.02+0,34*0,2=1,088

3).x3=0,6;   x2+1/2=0,5
    f12=z2=1,088                                        f22=2*0,4-
1,4084+1,088=0,4796
    y2+1/2=1,4084+1,088*0,1=1,5172        z2+1/2=1,088+0,4796*0,1=1,13596
    ?2=z2+1/2=1,13596                                ?2=2*0,5-
1,5172+1,13596=0,61876
    y3=1,4084+1,136*0,2=1,635592         z3=1,088+0,61876*0,2=1,211752

4).x4=0,8;   x3+1/2=0,7
    f13=z3=1,211752                                   f23=2*0,6-
1,636+1,212=0,77616
    y3+1/2=1,636+1,212*0,1=1,7567672     z3+1/2=1,212+0,776*0,1=1,289368
    ?3=z3+1/2=1,289368                               ?3=2*0,7-
1,7568+1,289=0,9326008
    y4=1,6+1,289*0,2=1,8934656              z4=1,212+0,93*0,2=1,39827216
5).x5=1;   y4+1/2=0,9
    f14=z4=1,39827216                                f24=2*0,8-
1,893+1,398=1,10480656
    y4+1/2=1,893+1,398*0,1=2,0332928
z4+1/2=1,398+1,105*0,1=1,508752816
    ?4=z4+1/2=1,508752816                          ?4=2*0,9-
2,03+1,5=1,27546
    y5=1,893+1,5*0,2=2,195216163           z5=1,398+1,275*0,2=1,65336416



3. Чтобы решить уравнение третьего порядка
                              y///=2x-y-y/+y//
на отрезке [0,1], с шагом h=0,2 и начальными условиями
                                   y0//=1
                                    y0/=1
                                    y0=1
необходимо сделать 3 замены:     y/=a                             y0/=a0=1
                                                         y//=a/=b
             y0//=b0=1
                                                         b/=2x-y-a+b

1).x1=0,2;   x1/2=0,1
      y(a1)=y(a0)+a0h
y(a1/2)=y(a0)+f10h/2
      a(b1)=a(b0)+?0h
a(b1/2)=a(b0)+f20h/2
      b(x1,y1,a1)=b(x0,y0,a0)+?0h
b(x1/2,y1/2,a1/2)=b(x0,y0,a0)+f30h/2
    f10=f(a0,y(a0))=1                                     y1/2=1+1*0,1=1,1
    f20=f(b0,a(b0))=1                                     a1/2=1+1*0,1=1,1
    f30=f(x0,y0,a0,b0)=-1                                b1/2=1-1*0,1=0,9
    ?0=a1/2=1,1
y(a1)=1+1,1*0,2=1,22
    ?0=b1/2=0,9
a(b1)=1+0,9*0,2=1,18
    ?0=2*0,1-1,1-1,1+0,9=-1,1                     b(x1,y1,a1)=1-
1,1*0,2=0,78

2).x2=0,4;   x1+1/2=x1+h/2=0,3
    f11=a1=1,18
y1+1/2=1,22+1,18*0,1=1.338
    f21=b1=0,78
a1+1/2=1,18+0,78*0,1=1,258
    f31=2*0,2-1,22-1,18+0,78=-1,22             b1+1/2=-1,22*0,1+0,78=0,658
    ?1=a1+1/2=1,258
y2=1,22+1,258*0,2=1,4716
    ?1=b1+1/2=0,658
a2=1,18+0,658*0,2=1,3116
    ?1=2*0,3-1,338-1,258+0,658=-1,338      b2=0,78-1,338*0,2=0,5124

3).x3=0,6;   x2+1/2=0,5
    f12=a2=1,3116
y2+1/2=1,47+1,3*0,1=1,60276
    f22=b2=0,5124
a2+1/2=1,3116+0,5*0,1=1.36284
    f32=2*0,4-1,47-1,31+0,512=-1,4708        b2+1/2=0,4-1,4*0,1=0,36542
    ?2=1,36284
y3=1,4716+1,3116*0,2=1,744168
    ?2=0,36542
a3=1,3116+0,3654*0,2=1,384664
    ?2=2*0,5-1,6-1,36+0,365=-1,60018         b3= 0,51-1,60018*0,2=0,192364

4).x4=0,8;   x3+1/2=0,7
    f13=1,384664
y3+1/2=1,74+1,38*0,1=1,8826364
    f23=0,192364
a3+1/2=1,38+0,19*0,1=1,4039204
    f33=2*0,6-1,7-1,38+0,19=-1,736488         b3+1/2=0,19-1,7*0,1=0,0187152


    ?3=1,4039204
y4=1,74+1,4*0,2=2,0249477
    ?3=0,0187152
a4=1,38+0,9187*0,2=1,388403
    ?3=2*0,7-1,88-1,4+0,0187=-1,8678416     b4=0,192-1,87*0,2=-0,1812235

5).x4=1;   x4+1/2=0,9
    f14=1,388403
y4+1/2=2,02+1,388*0,1=2,16379478
      f24=-0,1812235                                           a4+1/2=1,4-
0.181*0,1=1,370306608
    f34=2*0,8-2,02-1,388-0,18=-1,9945834      b4+1/2=-0,18-1,99*0,1=-
0,38066266
    ?4=1,3703
y5=2,02+1,37*0,2=2,2990038
    ?4=-0,38066                                                 a5=1,388-
0,38*0,2=1,3122669
    ?4=2*0,9-2,16-1,37-0,38=-2,114764056     b5=-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 – обращение к метке, задержка для просмотра результатов, очистка
        экрана, конец программы.


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