Метод милна для системы ОДУ |
1. Заголовок темы должен быть информативным. В противном случае тема удаляется ...
2. Все тексты программ должны помещаться в теги [code=pas] ... [/code], либо быть опубликованы на нашем PasteBin в режиме вечного хранения.
3. Прежде чем задавать вопрос, см. "FAQ", если там не нашли ответа, воспользуйтесь ПОИСКОМ, возможно такую задачу уже решали!
4. Не предлагайте свои решения на других языках, кроме Паскаля (исключение - только с согласия модератора).
5. НЕ используйте форум для личного общения, все что не относится к обсуждению темы - на PM!
6. Одна тема - один вопрос (задача)
7. Проверяйте программы перед тем, как разместить их на форуме!!!
8. Спрашивайте и отвечайте четко и по существу!!!
Метод милна для системы ОДУ |
Ozzя |
Сообщение
#1
|
Гуру Группа: Пользователи Сообщений: 1 220 Пол: Мужской Репутация: 16 |
ни у кого не завалялся?
|
Ozzя |
Сообщение
#2
|
Гуру Группа: Пользователи Сообщений: 1 220 Пол: Мужской Репутация: 16 |
Почитал. Если правильно понял, то систему ОДУ можно решать методом Рунге-Кутта, находя yi=f(xi) по методу Милна?
|
Sek |
Сообщение
#3
|
Гость |
Было бы не плохо отыскать бы программу на математике решения модели Мариока-Шимицу многошаговым методом Милна
|
Гость |
Сообщение
#4
|
Гость |
Program Miln;
Uses Crt; var x0,y0,z0,h:real; k1,k2,k3,k4:real; r1,r2,r3,r4:real; eps,abs_pogr:real; z,zkor,zpr,ypr,ykor,x,y:array [0..10] of real; i:integer; function f1(xa,ya,za:real):real; begin f1:=2*sqr(xa)+2*ya+za; end; function f2(xa,ya,za:real):real; begin f2:=1-2*sqr(xa)+2*ya-za; end; begin Clrscr; eps:=10e-6; x0:=0; y0:=1; z0:=1; h:=0.1; x[0]:=x0; y[0]:=y0; z[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],z[i]); r1:=h*f2(x[i],y[i],z[i]); k2:=h*f1(x[i]+h/2,y[i]+k1/2,z[i]+r1/2); r2:=h*f2(x[i]+h/2,y[i]+k1/2,z[i]+r1/2); k3:=h*f1(x[i]+h/2,y[i]+k2/2,z[i]+r2/2); r3:=h*f2(x[i]+h/2,y[i]+k2/2,z[i]+r2/2); k4:=h*f1(x[i]+h,y[i]+k3,z[i]+r3); r4:=h*f2(x[i]+h,y[i]+k3,z[i]+r3); y[i+1]:=y[i]+(k1+2*k2+2*k3+k4)/6; z[i+1]:=z[i]+(r1+2*r2+2*r3+r4)/6; x[i+1]:=x[i]+h; i:=i+1; end; i:=4; while x[i]<=1.0+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],z[i-3])-f1(x[i-2],y[i-2],z[i-2])+2*f1(x[i-1],y[i-1],z[i-1])); ykor[i]:=y[i-2]+(h/3)*(f1(x[i-2],y[i-2],z[i-2])+4*f1(x[i-1],y[i-1],z[i-1])+f1(x[i],ypr[i],z[i])); zpr[i]:=z[i-4]+(4*h)/3*(2*f2(x[i-3],y[i-3],z[i-3])-f2(x[i-2],y[i-2],z[i-2])+2*f2(x[i-1],y[i-1],z[i-1])); zkor[i]:=z[i-2]+(h/3)*(f2(x[i-2],y[i-2],z[i-2])+4*f2(x[i-1],y[i-1],z[i-1])+f2(x[i],zpr[i],z[i])); abs_pogr:=abs(ykor[i]-ypr[i])/29; if abs_pogr>eps then y[i]:=ykor[i] else y[i]:=ypr[i]; abs_pogr:=abs(zkor[i]-zpr[i])/29; if abs_pogr>eps then z[i]:=zkor[i] else z[i]:=zpr[i]; x[i+1]:=x[i]+h; i:=i+1; end; WriteLn; for i:=0 to 10 do begin WriteLn(x[i]:10:4,' ',y[i]:10:4,' ',z[i]:10:4); end; Readln; end. |
Текстовая версия | 19.03.2024 16:03 |