Прямые методы Рунге — Кутты

Семейство прямых методов Рунге — Кутты является обобщением метода Рунге — Кутты 4 порядка. Оно задаётся формулами

где — величина шага сетки по и вычисление нового значения проходит в этапов:

Конкретный метод определяется числом и коэффициентами и . Эти коэффициенты часто упорядочивают в таблицу (называемую таблицей Бутчера)

Для коэффициентов метода Рунге — Кутты должны быть выполнены условия для . Если требуется, чтобы метод имел порядок , то следует так же обеспечить условие

где — приближение, полученное по методу Рунге — Кутты. После многократного дифференцирования это условие преобразуется в систему полиномиальных уравнений на коэффициенты метода.

Листинг программы в приложении

Результат работы программы

С шагом 0, 2:

Vvedite znachenija koncov otrezka [a, b]

1 2

vveditenachalnoe znachenie funkcii y0 pri x=x0

Vvedite chislo znacheniy funkcii na promegutke (a, b]

Metod Eilera Metoda E-koshi Metod Runge-Koshi

x= 1.00 y1= 0.000000 y2= 0.000000 y3= 0.000000

x= 1.20 y1= 0.200000 y2= 0.216667 y3= 0.218779

x= 1.40 y1= 0.433333 y2= 0.467063 y3= 0.471049

x= 1.60 y1= 0.695238 y2= 0.746287 y3= 0.751990

x= 1.80 y1= 0.982143 y2= 1.050684 y3= 1.057997

x= 2.00 y1= 1.291270 y2= 1.377426 y3= 1.386272

Dlja prodolgenija nagmite lubuyu klavishu

С шагом 0, 1:

Vvedite znachenija koncov otrezka [a, b]

1 2

vveditenachalnoe znachenie funkcii y0 pri x=x0

Vvedite chislo znacheniy funkcii na promegutke (a, b]

Metod Eilera Metoda E-koshi Metod Runge-Koshi

x= 1.00 y1= 0.000000 y2= 0.000000 y3= 0.000000

x= 1.10 y1= 0.100000 y2= 0.104545 y3= 0.104841

x= 1.20 y1= 0.209091 y2= 0.218216 y3= 0.218785

x= 1.30 y1= 0.326515 y2= 0.340247 y3= 0.341073

x= 1.40 y1= 0.451632 y2= 0.469991 y3= 0.471060

x= 1.50 y1= 0.583891 y2= 0.606896 y3= 0.608197

x= 1.60 y1= 0.722817 y2= 0.750480 y3= 0.752005

x= 1.70 y1= 0.867993 y2= 0.900326 y3= 0.902067

x= 1.80 y1= 1.019052 y2= 1.056065 y3= 1.058015

x= 1.90 y1= 1.175666 y2= 1.217366 y3= 1.219521

x= 2.00 y1= 1.337543 y2= 1.383938 y3= 1.386293

Dlja prodolgenija nagmite lubuyu klavishu

Погрешность по методу Рунге < 0.000001

Погрешность по методу Рунге < 0.000002

Приложение

Листинг программы к заданию 1

program GausMethod;

uses Crt;

label

finish;

const

l=9;

type

matriceA=array[1..l, 1..l+1] of real;

vectorX=array[1..l] of real;

var

i, j, k, m, n: integer;

aa, bb, S: real;

a: matriceA;

x: vectorX;

ch: char;

{Podprorammy}

procedure PAUSA;

begin

writeln;

writeln('Dlja prodolgenija nagmite lubuju klavishu');

repeat

ch: =readkey

until (ch< > '');

end;

procedure OutputOfMatrice;

begin

ClrScr;

for i: = 1 to n do

begin

for j: = 1 to n+1 do

write (a[i, j]: 16: 7);

writeln;

end;

pausa;

end;

{Main Program}

begin

ClrScr;

writeln ('Vvedite chislo uravneniy n');

read(n);

writeln ('Vvedite rasshirennuju matricu');

writeln ('po odnomu elementu po strokam');

for i: = 1 to n do

begin

for j: = 1 to n+1 do

read(a[i, j]);

writeln;

end;

OutputOfMatrice;

for k: = 1 to n do

begin

{Poisk maksimalnogo elementa v k-toy stroke}

aa: =abs(a[k, k]);

i: =k;

for m: = k+1 to n do

begin

if abs(a[m, k])> aa then

begin

i: =m;

aa: =abs(a[m, k]);

end;

end;

if aa< 0.000001 then

begin

writeln;

writeln('Sistema ne imeet edinstvennogo reshenija');

pausa;

goto finish;

end;

{Perestanovka i-toy stroki,

sodergachey glavniy element s k-toy

strokoy}

if i< > k then

begin

for j: = k to n+1 do

begin

bb: =a[k, j];

a[k, j]: =a[i, j];

a[i, j]: =bb;

end;

end;

{Preobrazovanie k-toy stroki}

aa: =a[k, k];

a[k, k]: =1;

for j: = k+1 to n+1 do

a[k, j]: =a[k, j]/aa;

{Preobrazovanie s pomochju k-toj stroki ostalnih strok }

for i: = k+1 to n do

begin

bb: =a[i, k];

a[i, k]: =0;

if bb< > 0 then

begin

for j: = k+1 to n+1 do

a[i, j]: = a[i, j]-bb*a[k, j];

end;

end;

end;

OutputOfMatrice;

{koec prjamogo hoda v metode Gaussa}

for i: = n downto 1 do

begin

s: =a[i, n+1];

for j: = n downto i+1 do

s: =s-a[i, j]*x[j];

x[i]: =s;

end;

ClrScr;

for i: = 1 to n do

writeln('x', i, '=', x[i]);

pausa;

finish:;

end.

Листинг программы к заданию 2, 1

program metodPolovinngoDelenija;

uses crt;

label

finish;

var

a, b, d, e, x1, fa, f1: real;

ch: char;

{procedures/functions}

function f(x: real): real;

begin

f: =x*x*x-2*x*x-4*x+5;

end;

procedure OutputOfResult;

begin

clrscr;

writeln('x=', x1);

writeln('Pogreshnost rezultat =', d);

writeln ('Dlja prodolgenija nagmite lubuyu klavishu');

repeat ch: =readkey until ch< > '';

end;

{Main program}

begin

clrscr;

writeln('Vvedite po porjadku 3 chisla');

writeln('koordinaty koncov otrezka a, b, tochost epsilon');

readln(a, b, e);

if f(a)*f(b)> 0

then

goto finish;

repeat

x1: =(a+b)/2;

f1: =f(x1);

if f1=0 then

begin

d: =0;

OutputOfResult;

goto finish;

end;

fa: =f(a);

if fa*f1< 0 then

b: =x1

else

a: =x1;

d: =b-a;

OutputOfResult;

until d< e;

finish:

end.

Листинг программы к заданию 2, 2

program ReshenieUravnenijaMetodomIteraciy;

uses crt;

var

d, e, x, x1: real;

ch: char;

{procedures/functions}

function g(x: real): real;

begin

g: =(x*x*x-2*x*x+5)/4;

end;

procedure OutputOfResult;

begin

{clrscr; }

writeln('x=', x1);

writeln('Pogreshnost pribligenija d=', d);

writeln ('Dlja prodolgenija nagmite lubuyu klavishu');

repeat ch: =readkey until ch< > '';

end;

{Main program}

begin

clrscr;

writeln('Vvedite po porjadku 2 chisla');

writeln(' Nachalnoe znachenie i tochost epsilon/10');

readln(x, e);

repeat

x1: =g(x);

d: =abs(x1-x);

outputOfresult;

x: =x1;

until d< e;

end.

Листинг программы к заданию 3

uses crt;

const

l=15;

type

parameters=array [0..6] of real;

data=array [1..l] of real;

var

i, j, k, m: integer;

d1: real;

a, b, d: parameters;

x, y, x1, y1: data;

ch: char;

{Procedures}

procedure PAUSA;

begin

writeln;

writeln ('Dlja prodolgenija nagmite lubuyu klavishu');

repeat ch: =readkey until ch < > '';

end;

{Vichislenije parametrov lineinoi funkcii

pri viravnivanii eksperimentalnih dannyh}

procedure CalculOfParameters;

var

k0, b0, a1, b1, a2, b2, f2: real;

begin

a1: =0;

b1: =0;

a2: =0;

b2: =0;

for i: =1 to m do

begin

a1: =a1+x1[i];

b1: =b1+y1[i];

a2: =a2+x1[i]*x1[i];

b2: =b2+x1[i]*y1[i];

end;

d1: =m*a2-a1*a1;

k0: =(m*b2-a1*b1)/d1;

b0: =(b1*a2-a1*b2)/d1;

d1: =0;

f2: =0;

for i: =1 to m do

begin

f2: =f2+y1[i]*y1[i];

d1: =d1+(y1[i]-k0*x1[i]-b0)*(y1[i]-k0*x1[i]-b0);

end;

d[j]: =sqrt(d1/f2);

a[j]: =k0;

b[j]: =b0;

end;

{Main Program}

begin

clrscr;

writeln ('Vvedite chislo eksperimentalnih tochek m> 1');

read (m);

writeln ('Vvedite pari nacheniy xi, yi');

for i: =1 to m do

begin

read (x[i], y[i]);

{writeln(x[i], y[i]); }

end;

for j: =0 to 6 do

begin

case j of

0:

for i: =1 to m do

begin

x1[i]: =x[i];

y1[i]: =y[i];

end;

1:

for i: =1 to m do

y1[i]: =x[i]*y[i];

2:

for i: =1 to m do

y1[i]: =1/y[i];

3:

for i: =1 to m do

y1[i]: =x[i]/y[i];

4:

for i: =1 to m do

y1[i]: =ln(y[i]);

5:

for i: =1 to m do

begin

x1[i]: =ln(x[i]);

y1[i]: =y[i];

end;

6:

for i: =1 to m do

y1[i]: =y[i];

end;

CalculOfParameters;

end;

d1: =d[0];

k: =0;

for i: =1 to 6 do

if d[i]< d1 then

begin

d1: =d[i];

k: =i;

end;

write ('empirichskaja formula N',

k: 1, ' d=', d1: 8: 5);

writeln (' k=', a[k]: 8: 5, ' b=', b[k]: 8: 5);

pausa;

end.

Листинг программы к заданию 4, 1

program Integration;

{vich integrala po formule Simpsona}

uses crt;

var

i, n: integer;

a, b, h, s1, s2, s3, x, x1, x2: real;

ch: char;

{procedures}

function f(x: real): real;

begin

f: =ln(x+1);

end;

procedure PAUSA;

begin

writeln;

writeln ('Dlja prodolgenija nagmite lubuyu klavishu');

repeat ch: =readkey until ch < > '';

end;

{Main program}

begin

clrscr;

writeln('vvedite znachenija koncov otrezka [a, b]');

read (a, b);

writeln ('Vvedite chislo razbieniy otrezka n');

read(n);

h: =(b-a)/n;

x1: =a+h/2;

x2: =a;

s1: =f(x1);

s2: =(f(x2)+f(b))/2;

for i: =1 to n-1 do

begin

x1: =x1+h;

x2: =x2+h;

s1: =s1+f(x1);

s2: =s2+f(x2);

end;

s1: =h*s1;

s2: =h*s2;

s3: =(2*s1+s2)/3;

writeln('s1=', s1: 9: 6, ' s2=', s2: 9: 6, ' s3=', s3: 9: 6);

pausa;

end.

Листинг программы для задания 4, 2

program IntegrationbyGauss;

{vich integrala po metodu Gaussa}

uses crt;

var

i, n: integer;

a, b, c, c1, d, e, f1, f3, h, h1, s, s1, s2, x1, x2, x3: real;

ch: char;

{procedures}

function f(x: real): real;

begin

f: =ln(x+1);

end;

procedure PAUSA;

begin

writeln;

writeln ('Dlja prodolgenija nagmite lubuyu klavishu');

repeat ch: =readkey until ch < > '';

end;

{Main program}

begin

clrscr;

writeln('vvedite znachenija koncov otrezka [a, b]');

read (a, b);

writeln ('Vvedite tochnost` vichisleniya epsilon');

read(e);

c: =sqrt(3/5);

h1: =(b-a)/2;

c1: =c*h1;

x2: =(a+b)/2;

f1: =f(x2-c1);

f3: =f(x2+c1);

s1: =h1*(5*f1+8*f(x2)+5*f3)/9;

n: =2;

repeat

h: =(b-a)/n;

h1: =h/2;

c1: =c*h1;

x2: =a+h1;

X1: =X2-C1;

x3: =x2+c1;

s2: =0;

for i: =1 to n do

begin

s2: =s2+5*f(x1)+8*f(x2)+5*f(x3);

x1: =x1+h;

x2: =x2+h;

x3: =x3+h;

end;

s2: =s2*h1/9;

d: =abs(s1-s2)/63;

s1: =s2;

n: =2*n;

until d< e;

writeln ('Velichina integrala s=', s2: 9: 6,

' Pogreshnost` d=', d: 9: 6);

pausa;

end.

Листинг программы для задания 5

program DifEquationsOfirstOrder;

uses crt;

const

c: array[1..4] of real=(0, 0.5, 0.5, 1);

type

coef=array[0..4] of real;

var

i, j, m: integer;

a, b, h, x, y, y1, y2, y3: real;

k0, k: coef;

ch: char;

{Procedures}

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

begin

f: =(1-y+ln(x))/x;

end;

procedure PAUSA;

begin

writeln;

writeln ('Dlja prodolgenija nagmite lubuyu klavishu');

repeat ch: =readkey until ch < > '';

end;

{Main program}

begin

clrscr;

writeln('Vvedite znachenija koncov otrezka [a, b]');

read (a, b);

writeln('vveditenachalnoe znachenie funkcii y0 pri x=x0');

read (y);

writeln ('Vvedite chislo znacheniy funkcii na promegutke (a, b]');

read(m);

x: =a;

h: =(b-a)/m;

y1: =y;

y2: =y;

y3: =y;

writeln(' Metod Eilera Metoda E-koshi Metod Runge-Koshi');

writeln ('x=', x: 5: 2, ' y1=', y1: 9: 6, ' y2=', y2: 9: 6, ' y3=', y3: 9: 6);

for i: =1 to m do

begin

y1: =y1+h*f(x, y1); {Metod Eilera}

for j: =1 to 2 do

k0[j]: =h*f(x+2*c[j]*h, y2+2*c[j]*k0[j-1]);

y2: =y2+(k0[1]+k0[2])/2; {Metod Eilera-Koshi}

for j: =1 to 4 do {Metod runge-Kutta}

k[j]: =h*f(x+c[j]*h, y3+c[j]*k[j-1]);

y3: =y3+(k[1]+2*k[2]+2*k[3]+k[4])/6;

x: =x+h;

writeln ('x=', x: 5: 2, ' y1=', y1: 9: 6, ' y2=', y2: 9: 6, ' y3=', y3: 9: 6);

end;

pausa

end.