В этом методе начальные построения те же, только при вычислении площади криволинейной трапеции ординаты сверху соединяются ломаной линией.

Получается множество прямоугольных трапеций. Площадь одной трапеции равна:

Sтр=.h

Отсюда: y.h + .h + … +.h =

= h . + f(a + h) +…+ f(в-h)=+

Точность Е h2

Метод Симпсона (парабол)

В этом методе отрезок [a,в] разбивается на 2n частей, длинной h=и ординаты сверху соединяются кривой второго порядка (3 соседних точки).

Расчетная формула имеет вид :

у(y0 + 4y1+ 2y2 + 4y3 + …+ 4y2n*1 + y2n)

y0= f(a), y1= f(a+2h), y2= f(a+2h)… y2n-1= f(в-h)

y2n= f(в).

Для упрощения расчётов введём переменную сi

ci= , тогда формула примет вид:

у(f(a)+f(в)+(3+ci))

PROGRAM PRYAMOUGOLNIK;

CONST a= 0.4 ; b= 1.2; n=100;

var x,y:real;

i: integer;

function f (x:real):real;

begin

f:= (sin(x*x+0.5)/cos(x*x+0.5))/(1+2*x*x);

end;

begin

y:=0; x:=a;

for i:=0 to n do

begin

y:=y+f(x);

x:=x+(b-a)/n;

end;

y:=y*(b-a)/n;

writeln('y=',y:8:5);

readln;

end.

ОТВЕТ:

y=0.28099

PROGRAM TRAPEZIA;

CONST a=0.4; b=0.8; n=100;

var x,y,h,s:real;

n: integer;

function f (x:real):real;

begin

f:= (sin(x*x+0.5)/cos(x*x+0.5))/(1+2*x*x);

end;

begin

h:=(b-a)/n;

x:=a+h;

s:=(f(a)+f(b))/2;

while x<=b-h do

begin

s:=s+f(x);

x:=x+h;

end;

y:=s*h;

writeln('y=',y:8:5);

readln;

end.

ОТВЕТ:

y=0.28173

PROGRAM SIMPSON;

LABEL 10;

CONST a= 0.4; b=0.8; n=100;

var x,y,h,s:real;

c,n: integer;

function f (x:real):real;

begin

f:= (sin(x*x+0.5)/cos(x*x+0.5))/(1+2*x*x);

end;

begin

h:=(b-a)/(2*n);

s:=f(a)+f(b);

c:=1;

x:=a+h;

10: s:=s+(3+c)*f(x);

x:=x+h;

c:=-c;

if x<= b-h then goto 10;

y:=s*h/3;

writeln('y=',y:8:5);

readln;

end.

ОТВЕТ:

y= 0.27919

Решение интеграла в ППП "Eureka"

y=integ((sin(x^2+0.5)/cos(x^2+0.5))/(1+2*x^2),x,0.4,0.8)

y=0.2823564

2.4 Методы решения дифференциальных уравнений

При использовании различных протекающих во времени процессах первым этапом является составление дифференциального уравнения, описывающего процесс, а вторым - поиск решения этого уравнения. Дифференциальным уравнением называется равенство, связывающее значение аргумента, неизвестной функции некоторых ее производных. Наивысший порядок входящих в уравнение производных называется порядком дифференциального уравнения.

Рассмотрим уравнение вида:

y=f(x,y) (1)

Уравннение (1) имеет бесконечное множество решений (рис. 1) - через каждую точку плоскости проходит интегральная кривая. Чтобы выделить одну кривую, нужно указать точку плоскости, через которую проходит кривая, т.е. указать так называемые начальные уравнения (значения x=x0 и y=y0) (2)

Метод Эйлера

Одним из методов решения дифференциального уравнения (1) с начальным условием (2) является метод Эйлера.

Будем рассматривать уравнение (1) на некотором отрезке [a,b]. Пусть отрезок поделен на n частей с шагом .

Обозначим X0=a, X1=X0+h, X2=X0+2h,…, Xn=X0+nh=b. Обозначим искомые y(X1),…y(Xn) через y1…yn.

Методика решения уравнения (1) с начальными условиями (2) связяна с разложением решения в ряд Тейлора в h-окрестности точки X0.

При отбрасывании членов ряда, содержащие производные второго и высшего порядков, получим

где f(X,y) - правая часть уравнения (1).

Таким образом,

.

При достаточно малой величине шага h метод Эйлера дает решения с большей точностью, т.к. погрешность близка к h2 (h<<1) на каждом шаге интегрирования.

Метод Рунге-Кутта

Недостатком метода Эйлера является змедление вычислений при выборе малой величины шага h, задающей точность решения.

Наиболее распространенным методом численного интегрирования дифференциальных уравнений служит метод Рунге-Кутта, обеспечивающий ускорение за счет большей точности вычислений на каждом шаге. Точность метода Рунге-Кутта оценивается величиной E?h2.

Уточнение достигается за счет специального подбора координат четурех точек, в которых вычисляется первая производная f(x,y). Вместо первой производной h•f(x,y) используемой в формуле Эйлера, вычисляется усредненная первая производная fi.

Формулы интегрирования по методу Рунге-Кутта имеют вид:

где

h=(xn - x0)/n i=0,1,2,…n

y'=(1-y2)cos(x)+0.6y

при х0=0; xn=1; у0=0; h=0.1

program eyler;

label 100;

const h=0.1; x0=0; xk=1; у0=0;

х0=а;

var h,y,x:real;

i: integer;

function f (x,y: real): real;

begin

f:= (1-y*y)*cos(x)+0.6*y;

end;

begin

y:=y0; x:=x0;

100: y:=y+h*f(x,y);

x:=x+h;

writeln(` x=',x:5:1,' y=',y:8:5);

if x<xk then goto 100;

readln;

end.

Страницы: 1, 2, 3, 4