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

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

ДАЛПУ.

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

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

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

із курсу «Математичне моделювання на ЕОМ».

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

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

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

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

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

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

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

начало.

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

h, h/2.

k:=0.

xk+½:=xk+h/2 yk+½:=yk+f (xk, yk) h/2.

?k:= f (xk+½, yk+½) 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 рівних частин 17-ї та одержимо послідовність х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.

хк хк½ xk+h=xk1 х.

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

yNk/=yk+f (xk, yk) h/2=yk+yk+½ Отримуємо точку Nk/. У цьому точці будуємо таку дотичну: y (xk+½)=f (xk+½, yk+½)=?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) у неповній середній точці y/k+½=f (xk+½, yk+½) визначають ук+1.

Для оцінки похибки у точці хк проводять обчислення кк з кроком h, потім із кроком 2h і беруть 1/3 різниці цих значений:

| ук*-у (хк)|=1/3(yk*-yk), де у (х)-точное рішення диференціального уравнения.

Отже, методом Эйлера можна вирішити рівняння будь-яких порядків. Наприклад, щоб вирішити рівняння другого порядку y//=f (y/, y, x) з початковими умовами 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] з кроком h=(1−0)/5=0,2 Початкові умови: у0=1; Користуючись рекурентными формулами (4), знаходимо: 1). x1=0,2; х½=0,1; y (x1)=y (x0)+?0h; y (x½)=y (x0)+f (x0,y0)h/2; f (x0,y0)=2(0−1=-1 y (x½)=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+½=x1+h/2=0,2+0,1=0,3 y (x1+½)=y (x1)+f (x1,y (x1))h/2 f (x1,y1)=2(0,2−0,86=-0,46 y (x1+½)=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+½=x2+h/2=0,4+0,1=0,5 f (x2,y2)=2*0,4−0,8172=-0,0172 y2+½=0,8172−0,0172*0,1=0,81 548.

?2=2*0,5−0,81 548=0,18 452 y3=0,8172+0,18 452*0,2=0,854 104.

4).x4=0,8; x3+½=x3+h/2=0,6+0,1=0,7 f (x3,y3)=2*0,6−0,854 104=0,345 896 y3+½=0,854 104+0,345 896*0,1=0,8 886 936.

?3=2*0,7−0,89=0,5 113 064 y4=0,854 104+0,5 113 064*0,2=0,95 636 528.

5).x5=1; x4+½=0,8+0,1=0,9 f (x4,y4)=2*0,8−0,956=0,64 363 472 y4+½=0,956+0,643*0,1=1,20 728 752;

?4=2*0,9−1,02=0,779 271 248 y5=0,956+0,7792*0,2=1,11 221 953.

2. Дано рівняння другого порядку: y//=2x-y+y/.

Знаходимо рішення тому ж відрізку [0,1] з кроком h=0,2;

Заміна: y/=z z/=2x-y+z.

Початкові умови: у0=1 z0=1.

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

y (z½)=y (z0)+f (z0,y0)h/2 z (x½, y½)=z (x0,y0)+f (x0,y0,z0)h/2 f (z0,y0)=f10=1 f (x0,y0,z0)=f20=2*0- 1+1=0 y½=1+1*0,1=1,1 z½=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+½=0,3 f11=z1=1,02 f21=2*0,2- 1,2+1,02=0,22 y1+½=1,2+1,02*0,1=1,1 z1+½=1,02+0,22*0,1=1,042.

?1=z1+½=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+½=0,5 f12=z2=1,088 f22=2*0,4- 1,4084+1,088=0,4796 y2+½=1,4084+1,088*0,1=1,5172 z2+½=1,088+0,4796*0,1=1,13 596.

?2=z2+½=1,13 596 ?2=2*0,5- 1,5172+1,13 596=0,61 876 y3=1,4084+1,136*0,2=1,635 592 z3=1,088+0,61 876*0,2=1,211 752.

4).x4=0,8; x3+½=0,7 f13=z3=1,211 752 f23=2*0,6- 1,636+1,212=0,77 616 y3+½=1,636+1,212*0,1=1,7 567 672 z3+½=1,212+0,776*0,1=1,289 368.

?3=z3+½=1,289 368 ?3=2*0,7- 1,7568+1,289=0,9 326 008 y4=1,6+1,289*0,2=1,8 934 656 z4=1,212+0,93*0,2=1,39 827 216 5).x5=1; y4+½=0,9 f14=z4=1,39 827 216 f24=2*0,8- 1,893+1,398=1,10 480 656 y4+½=1,893+1,398*0,1=2,332 928 z4+½=1,398+1,105*0,1=1,508 752 816.

?4=z4+½=1,508 752 816 ?4=2*0,9- 2,03+1,5=1,27 546 y5=1,893+1,5*0,2=2,195 216 163 z5=1,398+1,275*0,2=1,65 336 416.

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; x½=0,1 y (a1)=y (a0)+a0h y (a½)=y (a0)+f10h/2 a (b1)=a (b0)+?0h a (b½)=a (b0)+f20h/2 b (x1,y1,a1)=b (x0,y0,a0)+?0h b (x½, y½, a½)=b (x0,y0,a0)+f30h/2 f10=f (a0,y (a0))=1 y½=1+1*0,1=1,1 f20=f (b0,a (b0))=1 a½=1+1*0,1=1,1 f30=f (x0,y0,a0,b0)=-1 b½=1−1*0,1=0,9.

?0=a½=1,1 y (a1)=1+1,1*0,2=1,22.

?0=b½=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+½=x1+h/2=0,3 f11=a1=1,18 y1+½=1,22+1,18*0,1=1.338 f21=b1=0,78 a1+½=1,18+0,78*0,1=1,258 f31=2*0,2−1,22−1,18+0,78=-1,22 b1+½=-1,22*0,1+0,78=0,658.

?1=a1+½=1,258 y2=1,22+1,258*0,2=1,4716.

?1=b1+½=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+½=0,5 f12=a2=1,3116 y2+½=1,47+1,3*0,1=1,60 276 f22=b2=0,5124 a2+½=1,3116+0,5*0,1=1.36 284 f32=2*0,4−1,47−1,31+0,512=-1,4708 b2+½=0,4−1,4*0,1=0,36 542.

?2=1,36 284 y3=1,4716+1,3116*0,2=1,744 168.

?2=0,36 542 a3=1,3116+0,3654*0,2=1,384 664.

?2=2*0,5−1,6−1,36+0,365=-1,60 018 b3= 0,51−1,60 018*0,2=0,192 364.

4).x4=0,8; x3+½=0,7 f13=1,384 664 y3+½=1,74+1,38*0,1=1,8 826 364 f23=0,192 364 a3+½=1,38+0,19*0,1=1,4 039 204 f33=2*0,6−1,7−1,38+0,19=-1,736 488 b3+½=0,19−1,7*0,1=0,187 152.

?3=1,4 039 204 y4=1,74+1,4*0,2=2,249 477.

?3=0,187 152 a4=1,38+0,9187*0,2=1,388 403.

?3=2*0,7−1,88−1,4+0,0187=-1,8 678 416 b4=0,192−1,87*0,2=-0,1 812 235.

5).x4=1; x4+½=0,9 f14=1,388 403 y4+½=2,02+1,388*0,1=2,16 379 478 f24=-0,1 812 235 a4+½=1,4- 0.181*0,1=1,370 306 608 f34=2*0,8−2,02−1,388−0,18=-1,9 945 834 b4+½=-0,18−1,99*0,1=- 0,38 066 266.

?4=1,3703 y5=2,02+1,37*0,2=2,2 990 038.

?4=-0,38 066 a5=1,388- 0,38*0,2=1,3 122 669.

?4=2*0,9−2,16−1,37−0,38=-2,114 764 056 b5=-0,181−2,1*0,2=-0,6 041 734.

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

uses crt, pram, kurs1₁; 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 (10 000); clrscr; end.

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

Програма перебуває у файлі kursova1. pas, і має 2 модуля, у яких містяться заставки. Модулі перебувають у файлах pram. tpu і kurs1₁.tpu. Для запуску файла kursova1. pas в Turbo Pascal необхідно натиснути F9. З’явиться перша заставка, далі натиснути enter і введення всі необхідні початкові умови: порядок похідною, коефіцієнти при членах рада, відтинок і початкові значення у (х0). На екрані виводиться крок обчислення і таблиця з відповідями. Після натискання enter виводиться друга заставка, після чого ж ми повертаємося до тексту программы.

ОПИС ПРОГРАММЫ.

1 — введення даних, які у програмі 2 — використання мітки, очищення екрана, введення вимог, рішення диференціального рівняння залежно від введення початкових умов 3 — присвоєння початкових умов диференціального рівняння третього порядку 4 — висновок таблиці зі значеннями 5 — введення формул методу Эйлера для рівняння третього порядка.

6 — присвоєння початкових умов рішення диференціального рівняння другого порядку 7 — висновок таблиці для рівняння другого порядку 8 — формули методу Эйлера для рівняння другого порядка.

9 — початкові умови для диференціального рівняння першого порядку 10 — формули методу Эйлера на вирішення рівняння першого порядку 11 — висновок таблиці 12 — звернення до мітці, затримка для перегляду результатів, очищення екрана, кінець программы.