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

[Tanechka]

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

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

Задание
Ток в электрической цепи описывается дифференциальным уравнением (см. стр. 24), которое необходимо решить, найдя y(0.1), y(0.2) ... y(1) с помощью метода Рунге-Кутта второго порядка с коррекцией в средней точке (в дифференциальном уравнении k = 2). По найденным значениям y определить количество тепла, выделяющееся на единичном сопротивлении за единицу времени, по формуле с помощью метода Симпсона.



program rg;

const
H=0.1;
K=2;
B=1;
A=0;
N=10;
X=0;
Y=1;

var
I {P, V, T, CH, Pt}: real;
DY1, DY2, DY3, DY4, DY5, DY6, DY7, DY8, DY9, DY10: real;

X10, X11, X12, X13, X20, X21, X22, X23, X30, X31, X32, X33, X40, X41, X42, X43, X50, X51, X52, X53, X60, X61, X62, X63, X70, X71, X72, X73, X80, X81, X82, X83, X90, X91, X92, X93, X100, X101, X102, X103, X200: real;

Y10, Y11, Y12, Y13, Y20, Y21, Y22, Y23, Y30, Y31, Y32, Y33, Y40, Y41, Y42, Y43, Y50, Y51, Y52, Y53, Y60, Y61, Y62, Y63, Y70, Y71, Y72, Y73, Y80, Y81, Y82, Y83, Y90, Y91, Y92, Y93, Y100, Y101, Y102, Y103, Y200: real;

K10, K11, K12, K13, K20, K21, K22, K23, K30, K31, K32, K33, K40, K41, K42, K43,K50, K51, K52, K53, K60, K61, K62, K63, K70, K71, K72, K73, K80, K81, K82, K83, K90, K91, K92, K93, K100, K101, K102, K103, K200: real;

begin
//
X10:=X;
X11:=X10+0.5*H;
X12:=X10+0.5*H;
X13:=X10+H;
X20:= X10+H;

X21:=X20+0.5*H;
X22:=X20+0.5*H;
X23:=X20+H;
X30:= X20+H;

X31:=X30+0.5*H;
X32:=X30+0.5*H;
X33:=X30+H;
X40:= X30+H;

X41:=X40+0.5*H;
X42:=X40+0.5*H;
X43:=X40+H;
X50:= X40+H;

X51:=X50+0.5*H;
X52:=X50+0.5*H;
X53:=X50+H;
X60:= X50+H;

X61:=X60+0.5*H;
X62:=X60+0.5*H;
X63:=X60+H;
X70:= X60+H;

X71:=X70+0.5*H;
X72:=X70+0.5*H;
X73:=X70+H;
X80:= X70+H;

X81:=X80+0.5*H;
X82:=X80+0.5*H;
X83:=X80+H;
X90:= X80+H;

X91:=X90+0.5*H;
X92:=X90+0.5*H;
X93:=X90+H;
X100:= X90+H;

X101:=X100+0.5*H;
X102:=X100+0.5*H;
X103:=X100+H;
X200:= X100+H;
// Y
Y10:=Y;
Y11:=Y10+0.5*Y10;
Y12:=Y11+0.5*Y11;
Y13:=Y12+0.5*Y12;

Y20:=Y+DY1;
Y21:=Y20+0.5*Y20;
Y22:=Y21+0.5*Y21;
Y23:=Y22+0.5*Y22;

Y30:=Y+DY2;
Y31:=Y30+0.5*Y30;
Y32:=Y31+0.5*Y31;
Y33:=Y32+0.5*Y32;

Y40:=Y+DY3;
Y41:=Y40+0.5*Y40;
Y42:=Y41+0.5*Y41;
Y43:=Y42+0.5*Y42;

Y50:=Y+DY4;
Y51:=Y50+0.5*Y50;
Y52:=Y51+0.5*Y51;
Y53:=Y52+0.5*Y52;

Y60:=Y+DY5;
Y61:=Y60+0.5*Y60;
Y62:=Y61+0.5*Y61;
Y63:=Y62+0.5*Y62;

Y70:=Y+DY6;
Y71:=Y70+0.5*Y70;
Y72:=Y71+0.5*Y71;
Y73:=Y72+0.5*Y72;

Y80:=Y+DY7;
Y81:=Y80+0.5*Y80;
Y82:=Y81+0.5*Y81;
Y83:=Y82+0.5*Y82;

Y90:=Y+DY8;
Y91:=Y90+0.5*Y90;
Y92:=Y91+0.5*Y91;
Y93:=Y92+0.5*Y92;

Y100:=Y+DY9;
Y101:=Y100+0.5*Y100;
Y102:=Y101+0.5*Y101;
Y103:=Y102+0.5*Y102;
Y200:=Y100+DY10;

// Ki
K10:=H*(1-sin(2*X10+Y10)+Y10/(2+X10));

K11:=H*(X11+H/2);
K12:=H*(X12+H/2);
K13:=H*(X13+H/2);

K20:=H*(1-sin(2*X20+Y20)+Y20/(2+X20));
K21:=H*(X21+H/2);
K22:=H*(X22+H/2);
K23:=H*(X23+H/2);

K30:=H*(1-sin(2*X30+Y30)+Y30/(2+X30));
K31:=H*(X31+H/2);
K32:=H*(X32+H/2);
K33:=H*(X33+H/2);

K40:=H*(1-sin(2*X40+Y40)+Y40/(2+X40));
K41:=H*(X41+H/2);
K42:=H*(X42+H/2);
K43:=H*(X43+H/2);

K50:=H*(1-sin(2*X50+Y50)+Y50/(2+X50));
K51:=H*(X51+H/2);
K52:=H*(X52+H/2);
K53:=H*(X53+H/2);

K60:=H*(1-sin(2*X60+Y60)+Y60/(2+X60));
K61:=H*(X61+H/2);
K62:=H*(X62+H/2);
K63:=H*(X63+H/2);

K70:=H*(1-sin(2*X70+Y70)+Y70/(2+X70));
K71:=H*(X71+H/2);
K72:=H*(X72+H/2);
K73:=H*(X73+H/2);

K80:=H*(1-sin(2*X80+Y80)+Y80/(2+X80));
K81:=H*(X81+H/2);
K82:=H*(X82+H/2);
K83:=H*(X83+H/2);

K90:=H*(1-sin(2*X90+Y90)+Y90/(2+X90));
K91:=H*(X91+H/2);
K92:=H*(X92+H/2);
K93:=H*(X93+H/2);

K101:=H*(1-sin(2*X101+Y101)+Y101/(2+X101));

K101:=H*(X101+H/2);
K102:=H*(X102+H/2);
K103:=H*(X103+H/2);
K200:=H*(1-sin(2*X103+Y103)+Y103/(2+X103));

// DY

DY1:=(1/6)*(K10+2*K11+2*K12+K13);
DY2:=(1/6)*(K20+2*K21+2*K22+K23);
DY3:=(1/6)*(K30+2*K31+2*K32+K33);
DY4:=(1/6)*(K40+2*K41+2*K43+K43);
DY5:=(1/6)*(K50+2*K51+2*K53+K53);
DY6:=(1/6)*(K60+2*K61+2*K63+K63);
DY7:=(1/6)*(K70+2*K71+2*K73+K73);
DY8:=(1/6)*(K80+2*K81+2*K83+K83);
DY9:=(1/6)*(K90+2*K91+2*K93+K93);
DY10:=(1/6)*(K100+2*K101+2*K102+K103);

//
writeln (' :');
writeln ('');

writeln ('---|','------|---------------|--------------|---------------------|--------------------|');
writeln (' |', ' |',' X |', ' Y |', ' Ki ', '| |');
writeln ('---|','------|---------------|--------------|---------------------|--------------------|');
writeln (' |',' 0 | ',X10, ' | ',Y, ' |',K10, ' | |');
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 1 | ',X11,' | ',Y11, ' | ',K11, ' | |');
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' 1 |',' 2 | ',X12,' | ',Y12, ' | ',K12, ' |', DY1,' |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 3 | ',X13,' | ',Y13, ' | ',K13, ' | |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 4 | ',X20,' | ',Y20:0:3, ' | | |' );
writeln ('___|','______|_______________|______________|_____________________|____________________|');

writeln (' |',' 1 | ',X21,' | ',Y21, ' | ',K21, ' | |');
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 2 | ',X22,' | ',Y22, ' | ',K22, ' |', DY2,' |' );
writeln (' 2 |','------|---------------|--------------|---------------------| |');
writeln (' |',' 3 | ',X23,' | ',Y23, ' | ',K23, ' | |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 4 | ',X30,' | ',Y30:0:3, ' | | |' );
writeln ('___|','______|_______________|______________|_____________________|____________________|');

writeln (' |',' 1 | ',X31,' | ',Y31, ' | ',K31, ' | |');
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 2 | ',X32,' | ',Y32, ' | ',K32, ' |', DY3,' |' );
writeln (' 3 |','------|---------------|--------------|---------------------| |');
writeln (' |',' 3 | ',X33,' | ',Y33, ' | ',K33, ' | |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 4 | ',X40,' | ',Y40:0:3, ' | | |' );
writeln ('___|','______|_______________|______________|_____________________|____________________|');

writeln (' |',' 1 | ',X41,' | ',Y41, ' | ',K41, ' | |');
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 2 | ',X42,' | ',Y42, ' | ',K42, ' |', DY4,' |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' 4 |',' 3 | ',X43,' | ',Y43, ' | ',K43, ' | |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 4 | ',X50,' | ',Y50:0:3, ' | | |' );
writeln ('___|','______|_______________|______________|_____________________|____________________|');

writeln (' |',' 1 | ',X51,' | ',Y51, ' | ',K51, ' | |');
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 2 | ',X52,' | ',Y52, ' | ',K52, ' |', DY5,' |' );
writeln (' 5 |','------|---------------|--------------|---------------------| |');
writeln (' |',' 3 | ',X53,' | ',Y53, ' | ',K53, ' | |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 4 | ',X60,' | ',Y60:0:3, ' | | |' );
writeln ('___|','______|_______________|______________|_____________________|____________________|');

writeln (' |',' 1 | ',X61,' | ',Y61, ' | ',K61, ' | |');
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 2 | ',X62,' | ',Y62, ' | ',K62, ' |', DY6,' |' );
writeln (' 6 |','------|---------------|--------------|---------------------| |');
writeln (' |',' 3 | ',X63,' | ',Y63, ' | ',K63, ' | |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 4 | ',X70,' | ',Y70:0:3, ' | | |' );
writeln ('___|','______|_______________|______________|_____________________|____________________|');

writeln (' |',' 1 | ',X71,' | ',Y71, ' | ',K71, ' | |');
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 2 | ',X72,' | ',Y72, ' | ',K72, ' |', DY7,' |' );
writeln (' 7 |','------|---------------|--------------|---------------------| |');
writeln (' |',' 3 | ',X73,' | ',Y73, ' | ',K73, ' | |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 4 | ',X80,' | ',Y80:0:3, ' | | |' );
writeln ('___|','______|_______________|______________|_____________________|____________________|');

writeln (' |',' 1 | ',X81,' | ',Y81, ' | ',K81, ' | |');
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 2 | ',X82,' | ',Y82, ' | ',K82, ' |', DY8,' |' );
writeln (' 8 |','------|---------------|--------------|---------------------| |');
writeln (' |',' 3 | ',X83,' | ',Y83, ' | ',K83, ' | |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 4 | ',X90,' | ',Y90:0:3, ' | | |' );
writeln ('___|','______|_______________|______________|_____________________|____________________|');
writeln (' |',' 1 | ',X91,' | ',Y91, ' | ',K91, ' | |');
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 2 | ',X92,' | ',Y92, ' | ',K92, ' |', DY9,' |' );
writeln (' 9 |','------|---------------|--------------|---------------------| |');
writeln (' |',' 3 | ',X93,' | ',Y93, ' | ',K93, ' | |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 4 | ',X100,' | ',Y100:0:3, ' | | |' );
writeln ('___|','______|_______________|______________|_____________________|____________________|');

writeln (' |',' 1 | ',X101,' | ',Y101, ' | ',K101, ' | |');
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 2 | ',X102,' | ',Y102, ' | ',K102, ' |', DY10,' |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln ('10 |',' 3 | ',X103,' | ',Y103, ' | ',K103, ' | |' );
writeln (' |','------|---------------|--------------|---------------------| |');
writeln (' |',' 4 | ',X200,' | ',Y200:0:3, ' | | |' );
writeln ('___|','______|_______________|______________|_____________________|____________________|');

writeln ('');
writeln ('');
writeln (' , ....');

I:=((B-A)/(3*N)*(Y*Y+Y200*Y200+4*(Y20*Y20+ Y40*Y40+ Y60*Y60+ Y80*Y80+ Y100*Y100)+2*(Y30*Y30+ Y50*Y50+ Y70*Y70+ Y90*Y90)));
{P:=(B-A)/(3*N);
V:=Y*Y+Y200*Y200;
T:=4*(Y20*Y20+ Y40*Y40+ Y60*Y60+ Y80*Y80+ Y100*Y100);
Ch:=2*(Y30*Y30+ Y50*Y50+ Y70*Y70+ Y90*Y90);

writeln (P);
writeln (V);
writeln (T);
writeln (Ch);}
writeln (' ', I:0:4);

end.

 

[Vladimir_S]

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

Выложите, пожалуйста, задание (со стр. 24).

 

[Tanechka]

По заданию методом рунге-кутта второго порядка.

 

[Vladimir_S]

Дорогая Танечка, что ж Вы тут за ужасов понаписали? Методом Рунге-Кутты второго порядка, он же метод Эйлера, Ваша задачка решается в 10-15 строчек!
Но есть вопросы.
1. "С коррекцией в средней точке" - это как? Или имеется в виду с пересчетом?
2. Как связать физическую задачу с математической? Что есть Y и x с точки зрения физики? Ток и время? Или еще что-то?

 

[Tanechka]

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

Полное задание вот:
Пересчет в средней точке.

 

[Vladimir_S]

Ну вот, пожалуйста, полное решение Вашей задачки. Как видите, вовсе не нужно создавать листинг размером с обеденный стол - всё куда проще:

Const
 h=0.1;
 N=10;

Var
 I:Array[0..N] of Real;
 j:integer;
 t,Z,Sum:Real;

Function F(x,y:real):real;
begin
 F:=1-Sin(x*2+y)+y/(x+2);
end;

Begin
 I[0]:=1;
 t:=0;
 writeln('t= ',t:3:1,'   I= ',I[0]:6:3);
 for j:=1 to N do
  begin
   t:=t+h;
   Z:=I[j-1]+h*F(t-h,I[j-1]);
   I[j]:=I[j-1]+h*(F(t-h,I[j-1])+F(t,Z))/2;
   writeln('t= ',t:3:1,'   I= ',I[j]:6:3);
  end;
 writeln;
 Sum:=(Sqr(I[0])+Sqr(I[N]))/3*h;
 for j:=1 to (N div 2) do
  Sum:=Sum+4.0*h/3.0*Sqr(I[2*j-1]);
 for j:=2 to (N div 2) do
  Sum:=Sum+2.0*h/3.0*Sqr(I[2*j-2]);
 Writeln('Heat: ', Sum:0:5);
 Readln
End.

P.S. Между прочим, правильно "Рунге-Кутты", а не "Кутта". Потому что фамилия математика - Кутта, а не Кутт.

 

[Tanechka]

МЕтод Рунге-Кутты

Спасибо большое) я учту на будущее обязательно)

 

[platon]

Владимир, помогите пожалуйста и мне, понимаю, что нет предела студенческой лени=) ибо дедлайн подходит к концу :tehnari_ru_837:, буду исправляться.



Ток в электрической цепи описывается дифференциальным уравнением (см. стр. 24), которое необходимо решить, найдя y(0.1), y(0.2) ... y(1) с помощью метода Рунге-Кутта второго порядка с коррекцией по средней производной (в дифференциальном уравнении k = 6 По найденным значениям y определить количество тепла, выделяющееся на единичном сопротивлении за единицу времени, по формуле (той же, что у Tanechka) с помощью метода трапеций.


уравнение, о котором сказано на странице 24 такое-же)

 

[Vladimir_S]

нет предела студенческой лени

Ну не настолько же!
Значит, так: до строки Sum:=... берем из этой темы (если формулы одинаковы), а начиная с Sum:=... - отсюда.

 

[platon]

значит у меня должно получиться так: (только возникает вопрос, куда подставлять k=6)?

Const
h=0.1;
N=10;

Var
I:Array[0..N] of Real;
j:integer;
t,Z,Sum:Real;

Function F(x,y:real):real;
begin
F:= 1-Sin(x*2+y)+y/(x+2);
end;

Begin
I[0]:=1;
t:=0;
writeln('t= ',t:3:1,' I= ',I[0]:6:3);
for j:=1 to N do
begin
t:=t+h;
Z:=I[j-1]+h*F(t-h,I[j-1]);
I[j]:=I[j-1]+h*(F(t-h,I[j-1])+F(t,Z))/2;
writeln('t= ',t:3:1,' I= ',I[j]:6:3);
end;
writeln;
Sum:=(Sqr(I[0])+Sqr(I[N]))/2*h;
for j:=1 to N-1 do
Sum:=Sum+h*Sqr(I[j]);
Writeln('Heat: ', Sum:0:5);
Readln
End.


P.S. проверьте, пожалуйста досконально, у вас, как у специалиста, это займет пру минут. у меня ставки высоки. спасибо

 

[Vladimir_S]

только возникает вопрос, куда подставлять k=6

Platon, ну что за?.. Смотрим внимательно задание:

Ценой неимоверного напряжения физических, интеллектуальных, а также астральных сил узреваем в аргументе синуса тот самый злосчастный параметр k. И отсюда выводим радикально-кардинальную разницу в задачах: если Танечке предписано подставить k=2, то Вам надлежит подставить k=6. И, следовательно, перекатывая кусок программы, сделанной мною для Танечки, Вам нужно было всё-таки внести некоторые коррективы, а именно в Вашем случае функция будет такой:

Function F(x,y:real):real;
begin
F:= 1-Sin(x*6+y)+y/(x+2);
end;

Исправляем и радуемся жизни.

 

[Vladimir_S]

Victor81, для Вашей задачи создана отдельная тема: http://www.tehnari.ru/f43/t96383/ .

 

[platon]

да вопрос был действительно глупый, я когда переписывал сам додумал. спасибо большое)

 

[Vladimir_S]

да вопрос был действительно глупый, я когда переписывал сам додумал. спасибо большое)

Между прочим, если Вы из одной команды с Cbvrf86 (а я сильно подозреваю, что это так и есть), то, как сейчас обнаружилось, преподы требуют применить несколько иную модификацию метода Рунге-Кутты. Посмотрите эту тему: http://www.tehnari.ru/f43/t96332/.