double arrow

Второй усовершенствованный Эйлера.

Министерство образования Российской Федерации

ФГБОУ ВПО Ивановский государственный химико-технологический университет

Кафедра прикладной математики

 

 

ЛАБОРАТОРНАЯ РАБОТА № 5

По дисциплине: «Математические методы и модели в расчётах на ЭВМ»

Тема: Приближенное решение обыкновенных дифференциальных уравнений и систем обыкновенных дифференциальных уравнений.

Вариант 3

 

Выполнил: ст. группы 2.36

А.А. Машкова__________

Проверил: ст.преп.

С.В. Кулакова___________

 

 

Иваново 2011

Задание:

1)       Решить дифференциальное уравнение x^2*y’+y^2=0 на отрезке [-1;-0.1] с шагом h=0.1 методом Милна 4-го порядка и заданными начальными условиями y(-1)=1.

Решить систему ОДУ

  y’=z

  z’=-0.33xz-y

на отрезке [0;1] с шагами h1=0.05, h2=0.1 вторым усовершенствованным методом эйлера с заданными начальными условиями: y(0)=0, z(0)=0.25.

Теоретическая часть

Метод Милна 4-го порядка

 

2-ой усовершенствованный метод Эйлера (метод Рунге-Кутта 2-го порядка)

Рассмотрим метод Рунге-Кутта второго порядка. В этом методе величины y(i+1) вычисляются по следующим формулам:

Y(i+1) = y(i) + Δy(i)

 Δy(i)=Δy(i1)+Δy(i2)

 (см рис 6,6)

 

Тогда

Обозначим , тогда

Геометрически это означает, что мы определяем направление интегральной кривой в исходной точке (xi, yi) и во вспомогательной точке , а в качестве окончательного берем среднее из этих направлений.

Локальная погрешность третьего порядка:

e=у’’’*h^3/6

Глобальная погрешность второго порядка:O(h^2)

Оценка погрешности с двойным просчетом: (|I(h)-I(h/2)|)/3

 

Рассчетная часть

Метод Милна 4-го порядка

Определение начального отрезка производится по формуле Рунге-Кутта:

Сначала нужно преобразовать уравнение x^2*y’+y^2=0 в y’=-(y^2/x^2)

K1=h*(-(y^2/x^2))=0.1*(-(1)^2/(-1)^2)=-0.1

K2= h*(-((y+K1/2)^2)/((x+h/2)^2))=0.1*(-(1+(-0.1)/2)^2/((-1+0.1/2)^2))=-0.1

K3= h*(-((y+K2/2)^2)/((x+h/2)^2))= 0.1*(-(1+(-0.1)/2)^2/((-1+0.1/2)^2))=-0.1

K4= h*(-((y+K3)^2)/((x+h)^2))=0.1*(-((1-0.1)^2)/((-1+0.1)^2)=-0.1

Q=abs((K2-K3)/(K1-K2))=0

 

B3: =B2+D3

C3: =C2+(E2+2*H2+2*J2+L2)/6

E2: =-D2*((C2*C2)/(B2*B2))

F2: =B2+D2/2

G2: =C2+E2/2

H2: =-D2*((G2*G2/F2/F2))

I2: =C2+H2/2

J2: =-D2*((G2*G2/F2/F2))

K2: =C2+J2

L2: =-D2*((K2*K2/B3/B3))

M2: =ABS((H2-J2)/(E2-H2))

Остальные значения находим по методу Милна 4-го порядка:

 

 

Погрешность выбираем максимальную, т.е. e=0,01964

 

D10: =C2+4/3*(2*E3-E4+2*E5)

D11: =C3+4*(2*D4-D5+2*G10)/3

D12: =C4+4*(2*E5-G10+2*G11)/3

D13: =C5+4*(2*G10-G11+2*G12)/3

D14: =F10+4*(2*G11-G12+2*G13)/3

E10: =-C10*D10*D10/B10/B10

F10: =C4+(E4+4*E3+E10)/3

F11: =C5+(E5+4*G10+E11)/3

F12: =F10+(G10+4*G11+E12)/3

G10:= =-C10*F10*F10/B10/B10

H10: =ABS(F10-D10)/29

 

Программа на Паскале.

Program Miln;

Uses Crt; const n=10;

Var

  x0,y0,z0,h:real;

k1,k2,k3,k4:real;

eps,abs_pogr:real;

ypr,ykor,x,y: array [0..n] of real;

i:integer;

function f1(xa,ya:real):real;

Begin

  f1:=-ya*ya/xa/xa;

end;

Begin

  Clrscr;

eps:=0.001;

x0:=-1;

y0:=1;

h:=0.1;

x[0]:=x0;

y[0]:=y0;

i:=0;

{ Gotovim 1-e 3 tochki po metodu runge-kutta }

while i<=3 do

begin

      k1:=h*f1(x[i],y[i]);

k2:=h*f1(x[i]+h/2,y[i]+k1/2);

k3:=h*f1(x[i]+h/2,y[i]+k2/2);

k4:=h*f1(x[i]+h,y[i]+k3);

y[i+1]:=y[i]+(k1+2*k2+2*k3+k4)/6;

x[i+1]:=x0+i*h;

i:=i+1;

end;

i:=4;

while x[i]<=x0+i*h{eps> abs_pogr} do

begin

  { etap prognoza i korrektsii}

ypr[i]:=y[i-4]+(4*h)/3*(2*f1(x[i-3],y[i-3])-f1(x[i-2],y[i-2])+2*f1(x[i-1],y[i-1]));

 

ykor[i]:=y[i-2]+(h/3)*(f1(x[i-2],y[i-2])+4*f1(x[i-1],y[i-1])+f1(x[i],ypr[i]));

abs_pogr:=abs(ykor[i]-ypr[i])/29;

if abs_pogr>eps then

        y[i]:=ykor[i]

Else

          y[i]:=ypr[i];

     i:=i+1;

end;

WriteLn;

for i:=0 to 10 do

begin

      WriteLn(x[i]:10:4,' ',ykor[i]:10:4,' ');

end;

Readln;

end.

Второй усовершенствованный Эйлера.

Погрешность высчитывается с помощью двойного просчета:

Выбираем максимальную погрешность:e=

h=0.1

F0=z0=0.25

V0=-0.33*x0*z0-y0=-0.33*0*0.25-0=0

Y’1=y0+h*f0=0+0.1*0.25=0.025

Z’1=z0+h*v0=0.25+0.1*0=0.25

F’1=z’1=0.25

V’1=-0.33*x1*z’1-y’1=-0.03325

Y1=y0+h*(f+f’)/2=0.025

Z1=z0+h*(v+v’)/2=0.2483375

F1=z1=0.248338

V1=-0.33*x1*z1-y1=-0.0332

Y’2=y1+h*f1=0+0.1*0.25=0.049834

Z’2=z1+h*v1=0.25+0.1*0=0.245018

F’2=z’2=0.245018

V’2=-0.33*x2*z’2-y’2=-0.066

Y2=y0+h*(f+f’)/2=0.049668

Z2=z1+h*(v+v’)/2=0.2433775

Аналогичные вычисления проводятся с h=0.05

F0=0,25

V0=0

Y’1=0,0125

Z’1=0,25

F’1=0,25

V’1=-0,01663

Y1=0

Z1=0,25

F1=0,249584

V1=-0,01662

Y’2=0,024979

Z’2=0,248753

F’2=0,248753

V’2=-0,03319

Y2=0,0125

Z2=0,2495844

Вычисление погрешности: е=|Y(h)-Y(h/2)|/3

E1(0)=0-0=0

E2(0)=0.25-0.25=0

E1(1)=0.05-0.049668=0.000111

E2(1)=0.24335-0.2433775= 9,16542E-06

Из погрешностей выбираем максимальную для данной системы ОДУ::e=

Формулы в excel:

 

С3: =C2+E3*(F2+J2)/2

D3: =D2+E3*(G2+K2)/2

F2: =D2

G2: =-0,33*B2*D2-C2

H2: =C2+E2*F2

I2: =D2+E2*G2

J2: =I2

K2: =-0,33*B3*I2-H2

M2: =ABS(C14-C2)/3

N2: =ABS(D14-D2)/3

 

Программа на Паскале

program eiler;

function f(z:real):real;

Begin

  f:=z;

end;

function p(x:real;y:real;z:real):real;

Begin

  p:=-0/33*x*z-y;

end;

Var

  h:real;

 x:real;

 y1:real;

 y2:real;

 z1:real;

 z2:real;

Begin

  h:=0.1;

  while h<=0.2 do

 begin

  y1:=0;

z1:=0.25;

x:=0;

writeln('h= ',h:1:1);

while x<1 do

begin

   y2:=y1+h*f(z1);

z2:=z1+h*p(x,y1,z1);

y2:=y1+h/2*(f(z1)+f(z2));

z2:=z1+h/2*(p(x,y1,z1)+p(x,y2,z2));

z1:=z2;

y1:=y2;

x:=x+h;

writeln('x= ',x:2:4,' y= ',y2:2:4,' z= ',z2:2:4);

end;

h:=h+0.1;

  end;

 readln;

end.

 

 

Вывод:

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

 


Понравилась статья? Добавь ее в закладку (CTRL+D) и не забудь поделиться с друзьями:  



Сейчас читают про: