Скачать .docx | Скачать .pdf |
Реферат: Численный расчет дифференциальных уравнений
Міністерство освіти України
ДАЛПУ
Кафедра автоматизації
технологічних процесів і приладобудування
КУРСОВА РОБОТА
з курсу “Математичне моделювання на ЕОМ”
на тему “Розв’язок диференціального рівняння
виду ап у(п) +ап-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 +αk h
нет 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,…); каждое звено Мi Mi+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 )+α0 h; 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 )+α1 h; 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 )+α0 h z(x1 ,y1 )=z(x0 ,y0 )+β0 h
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 )+a0 h y(a1/2 )=y(a0 )+f10 h/2
a(b1 )=a(b0 )+β0 h a(b1/2 )=a(b0 )+f20 h/2
b(x1 ,y1 ,a1 )=b(x0 ,y0 ,a0 )+γ0 h b(x1/2 ,y1/2 ,a1/2 )=b(x0 ,y0 ,a0 )+f30 h/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 – обращение к метке, задержка для просмотра результатов, очистка
экрана, конец программы.