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

                          ДАЛПУ

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

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

            

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

      

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

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

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

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

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

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

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

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

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

                                                                                                                   

                                       

1998

                                                                       

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

                                            

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

начало

у^/=f(x,y)

y(x₀)=y₀

                                              x₀, x₀+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+α_kh

                                    _нет         k:=n

                                                                                    _да

 

                                                x₀, y₀,

x1, y1…

                                                x_n, y_n

                                                конец

           

             

               

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

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

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

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

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

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

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

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

                                     x=x₀, y(x₀)=y₀              (2)

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

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

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

При этом искомая интегральная кривая у=у(х), проходящая через точку М₀(х₀, у₀), заменяется ломаной М₀М₁М₂… с вершинами М_i(x_i, y_i) (i=0,1,2,…); каждое звено М_iM_i+1 этой ломаной, называемой ломаной Эйлера, имеет направление, совпадающее с направлением той интегральной кривой уравнения (1), которая проходит через точку М_i.

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

                           |f(x, y₁)- f(x, y₂)| £ N|y₁-y₂|  (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₀)=y₀. Разобьем наш участок интегрирования на n                                          

                                          равных частей. На малом участке [x₀,x₀+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₀)=y^/₀, y(x₀)=y₀, выполняется замена:

                                            y^/=z

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

Тем самым преобразуются начальные условия: y(x₀)=y₀, z(x₀)=z₀, z₀=y^/₀.

                                                                 
                    

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

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

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

y^/=2x-y

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

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

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

1). x₁=0,2;  х_1/2=0,1;     y(x₁)=y(x₀)+α₀h;   y(x_1/2)=y(x₀)+f(x₀,y₀)h/2;

      f(x₀,y₀)=2*0-1=-1  

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

      α₀=2*0,1-0,9=-0,7  

       y₁=1-0,1*0,2=0,86

2). y(x₂)=y(x₁)+α₁h;   x₂=0,2+0,2=0,4;   x_1+1/2=x₁+h/2=0,2+0,1=0,3

     y(x_1+1/2)=y(x₁)+f(x₁,y(x₁))h/2  

     f(x₁,y₁)=2*0,2-0,86=-0,46  

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

     α₁=2*0,3-0,814=-0,214

     y₂=0,86-0,214*0,2=0,8172

3). x₃=0,4+0,2=0,6;   x_2+1/2=x₂+h/2=0,4+0,1=0,5

     f(x₂,y₂)=2*0,4-0,8172=-0,0172 

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

     α₂=2*0,5-0,81548=0,18452 

     y₃=0,8172+0,18452*0,2=0,854104

 4).x₄=0,8;   x_3+1/2=x₃+h/2=0,6+0,1=0,7

     f(x₃,y₃)=2*0,6-0,854104=0,345896

     y_3+1/2=0,854104+0,345896*0,1=0,8886936

     α₃=2*0,7-0,89=0,5113064

     y₄=0,854104+0,5113064*0,2=0,95636528

 

5).x₅=1;   x_4+1/2=0,8+0,1=0,9

    f(x₄,y₄)=2*0,8-0,956=0,64363472

    y_4+1/2=0,956+0,643*0,1=1,020728752;  

    α₄=2*0,9-1,02=0,779271248

    y₅=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

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

                                           z₀=1

1).x₁=0,2;   x_1/2=0,1

   y(z₁)=y(z₀)+α₀h                                   z(x₁,y₁)=z(x₀,y₀)+β₀h              

   y(z_1/2)=y(z₀)+f(z₀,y₀)h/2                      z(x_1/2,y_1/2)=z(x₀,y₀)+f(x₀,y₀,z₀)h/2

   f(z₀,y₀)=f₁₀=1                                       f(x₀,y₀,z₀)=f₂₀=2*0-1+1=0

   y_1/2=1+1*0,1=1,1                                 z_1/2=1+0*0,1=1

   α₀=z₀=1                                                β₀=2*0,1-1,1+1=0,1

   y₁=1+0,2*1=1,2                                   z₁=1+0,2*0,1=1,02

 

2).x₂+0,4;   x_1+1/2=0,3

    f₁₁=z₁=1,02                                         f₂₁=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

    α₁=z_1+1/2=1,042                                   β₁=2*0,3-1,302+1,042=0,34

    y₂=1,2+1,042*0,2=1,4084                  z₂=1.02+0,34*0,2=1,088

3).x₃=0,6;   x_2+1/2=0,5

    f₁₂=z₂=1,088                                        f₂₂=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

    α₂=z_2+1/2=1,13596                                β₂=2*0,5-1,5172+1,13596=0,61876

    y₃=1,4084+1,136*0,2=1,635592         z₃=1,088+0,61876*0,2=1,211752

4).x₄=0,8;   x_3+1/2=0,7

    f₁₃=z₃=1,211752                                   f₂₃=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

    α₃=z_3+1/2=1,289368                               β₃=2*0,7-1,7568+1,289=0,9326008

    y₄=1,6+1,289*0,2=1,8934656              z₄=1,212+0,93*0,2=1,39827216

5).x₅=1;   y_4+1/2=0,9

    f₁₄=z₄=1,39827216                                f₂₄=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

    α₄=z_4+1/2=1,508752816                          β₄=2*0,9-2,03+1,5=1,27546

    y₅=1,893+1,5*0,2=2,195216163           z₅=1,398+1,275*0,2=1,65336416

   

   

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

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

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

y₀^//=1

y₀^/=1

y₀=1

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

                                                         y^//=a^/=b                        y₀^//=b₀=1

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

1).x₁=0,2;   x_1/2=0,1

      y(a₁)=y(a₀)+a₀h                                        y(a_1/2)=y(a₀)+f₁₀h/2

      a(b₁)=a(b₀)+β₀h                                        a(b_1/2)=a(b₀)+f₂₀h/2

      b(x₁,y₁,a₁)=b(x₀,y₀,a₀)+γ₀h                       b(x_1/2,y_1/2,a_1/2)=b(x₀,y₀,a₀)+f₃₀h/2

    f₁₀=f(a₀,y(a₀))=1                                     y_1/2=1+1*0,1=1,1

    f₂₀=f(b₀,a(b₀))=1                                     a_1/2=1+1*0,1=1,1

    f₃₀=f(x₀,y₀,a₀,b₀)=-1                                b_1/2=1-1*0,1=0,9

    α₀=a_1/2=1,1                                              y(a₁)=1+1,1*0,2=1,22

    β₀=b_1/2=0,9                                             a(b₁)=1+0,9*0,2=1,18

    γ₀=2*0,1-1,1-1,1+0,9=-1,1                     b(x₁,y₁,a₁)=1-1,1*0,2=0,78

2).x₂=0,4;   x_1+1/2=x₁+h/2=0,3

    f₁₁=a₁=1,18                                              y_1+1/2=1,22+1,18*0,1=1.338

    f₂₁=b₁=0,78                                              a_1+1/2=1,18+0,78*0,1=1,258

    f₃₁=2*0,2-1,22-1,18+0,78=-1,22             b_1+1/2=-1,22*0,1+0,78=0,658

    α₁=a_1+1/2=1,258                                        y₂=1,22+1,258*0,2=1,4716

    β₁=b_1+1/2=0,658                                        a₂=1,18+0,658*0,2=1,3116

    γ₁=2*0,3-1,338-1,258+0,658=-1,338      b₂=0,78-1,338*0,2=0,5124

3).x₃=0,6;   x_2+1/2=0,5

    f₁₂=a₂=1,3116                                           y_2+1/2=1,47+1,3*0,1=1,60276       

    f₂₂=b₂=0,5124                                           a_2+1/2=1,3116+0,5*0,1=1.36284

    f₃₂=2*0,4-1,47-1,31+0,512=-1,4708        b_2+1/2=0,4-1,4*0,1=0,36542

    α₂=1,36284                                               y₃=1,4716+1,3116*0,2=1,744168

    β₂=0,36542                                               a₃=1,3116+0,3654*0,2=1,384664

    γ₂=2*0,5-1,6-1,36+0,365=-1,60018         b₃= 0,51-1,60018*0,2=0,192364

4).x₄=0,8;   x_3+1/2=0,7

    f₁₃=1,384664                                             y_3+1/2=1,74+1,38*0,1=1,8826364

    f₂₃=0,192364                                             a_3+1/2=1,38+0,19*0,1=1,4039204

    f₃₃=2*0,6-1,7-1,38+0,19=-1,736488         b_3+1/2=0,19-1,7*0,1=0,0187152

   

    α₃=1,4039204                                             y₄=1,74+1,4*0,2=2,0249477

    β₃=0,0187152                                             a₄=1,38+0,9187*0,2=1,388403 

    γ₃=2*0,7-1,88-1,4+0,0187=-1,8678416     b₄=0,192-1,87*0,2=-0,1812235

5).x₄=1;   x_4+1/2=0,9

    f₁₄=1,388403                                               y_4+1/2=2,02+1,388*0,1=2,16379478

      f₂₄=-0,1812235                                           a_4+1/2=1,4-0.181*0,1=1,370306608

    f₃₄=2*0,8-2,02-1,388-0,18=-1,9945834      b_4+1/2=-0,18-1,99*0,1=-0,38066266

    α₄=1,3703                                                    y₅=2,02+1,37*0,2=2,2990038

    β₄=-0,38066                                                 a₅=1,388-0,38*0,2=1,3122669  

    γ₄=2*0,9-2,16-1,37-0,38=-2,114764056     b₅=-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 и ввести все необходимые начальные условия: порядок производной, коэффициенты при членах рада, отрезок и начальные значения у(х₀). На экране выводится шаг вычисления и таблица с ответами. После нажатия enter выводится вторая заставка, после чего мы возвращаемся к тексту программы. 

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

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

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

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

      условий

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

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

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

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

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

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

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

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

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

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

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

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

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

              Міністерство освіти України                           ДАЛПУ                Кафедра автоматизації   технологічних процесів і приладобудування                            КУРСОВА    РОБОТА                 з курсу “Матема

 

• 14:21 05 Июля 2016
  Российско-китайский медуниверситет откроет магистерскую программу по TCM
• 01:21 05 Июля 2016
  На форуме МГУ Россия и Китай обсудят научно-образовательное сотрудничество
• 11:00 04 Июля 2016
  Реморенко: "ЕГЭ-туризм" в России сходит на нет
• 10:52 04 Июля 2016
  Летняя российско-австрийская школа откроется в Москве
• 15:31 01 Июля 2016
  В СПбГУ могут воссоздать географический факультет

Заказать уникальную работу