Решение дифференциального уравнения

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

По теме: "Решение дифференциального уравнения"

1. Постановка задачи и метод решения

Решить дифференциальное уравнение у^/=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) линией. Получаем точку М_к(х_к, у_к).

Через М_к проводим касательную: у=у_к=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^/₀.

2. Решение контрольного примера

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

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₀)=20-1=-1

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

?₀=20,1-0,9=-0,7

y₁=1-0,10,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₁)=20,2-0,86=-0,46

y(x_1+1/2)=0,86-0,460,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.

3. Запуск программы на выполнение

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

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

4. Описание программы

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

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

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

условий

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

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

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

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

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

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

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

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

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

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

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

12 - обращение к метке, задержка для просмотра результатов, очистка экрана, конец программы.



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