MATLAB - Адаптивный размер шага Рунге-Кутта
Я запрограммировал в MATLAB адаптивный размер шага RK4 для решения системы ODE. Код выполняется без ошибок, однако он не дает желаемой кривой, когда я пытаюсь построить x против y. Вместо того, чтобы иметь форму тороида, я просто получаю плоскую линию. Это видно из того факта, что r выводит постоянное значение. После проверки выходных данных каждой строки они не выводят константы или ошибки или inf или NaN, а выводят как действительный, так и мнимый компонент (комплексные числа). Я понятия не имею, почему это происходит, и я считаю, что это является источником моей проблемы.
function AdaptRK4()
parsec = 3.08*10^18;
r_1 = 8.5*1000.0*parsec; % in cm
theta_1 = 0.0;
a = 0.5*r_1;
gam = 1;
grav = 6.6720*10^-8;
amsun = 1.989*10^33;
amg = 1.5d11*amsun;
gm = grav*amg;
u_1 = 20.0*10^5;
v = sqrt(gm/r_1);
time = 0.0;
epsilon = 0.00001;
m1 = 0.5;
m2 = 0.5;
m3 = 0.5;
i = 1;
nsteps = 50000;
deltat = 5.0*10^12;
angmom = r_1*v;
angmom2 = angmom^2.0;
e = -2*10^5.0*gm/r_1+u_1*u_1/2.0+angmom2/(2.0*r_1*r_1);
for i=1:nsteps
deltat = min(deltat,nsteps-time);
fk3_1 = deltat*u_1;
fk4_1 = deltat*(-gm*r_1*r_1^(-gam)/(a+r_1)^(3- gam)+angmom2/(r_1^3.0));
fk5_1 = deltat*(angmom/(r_1^2.0));
r_2 = r_1+fk3_1/4.0;
u_2 = u_1+fk4_1/4.0;
theta_2 = theta_1+fk5_1/4.0;
fk3_2 = deltat*u_2;
fk4_2 = deltat*(-gm*r_2*r_2^(-gam)/(a+r_2)^(3-gam)+angmom2/(r_2^3.0));
fk5_2 = deltat*(angmom/(r_2^2.0));
r_3 = r_1+(3/32)*fk3_1 + (9/32)*fk3_2;
u_3 = u_1+(3/32)*fk4_1 + (9/32)*fk4_2;
theta_3 = theta_1+(3/32)*fk5_1 + (9/32)*fk5_2;
fk3_3 = deltat*u_3;
fk4_3 = deltat*(-gm*r_3*r_3^(-gam)/(a+r_3)^(3-gam)+angmom2/(r_3^3.0));
fk5_3 = deltat*(angmom/(r_3^2.0));
r_4 = r_1+(1932/2197)*fk3_1 - (7200/2197)*fk3_2 + (7296/2197)*fk3_3;
u_4 = u_1+(1932/2197)*fk4_1 - (7200/2197)*fk4_2 + (7296/2197)*fk4_3;
theta_4 = theta_1+(1932/2197)*fk5_1 - (7200/2197)*fk5_2 + (7296/2197)*fk5_3;
fk3_4 = deltat*u_4;
fk4_4 = deltat*(-gm*r_4*r_4^(-gam)/(a+r_4)^(3-gam)+angmom2/(r_4^3.0));
fk5_4 = deltat*(angmom/(r_4^2.0));
r_5 = r_1+(439/216)*fk3_1 - 8*fk3_2 + (3680/513)*fk3_3 - (845/4104)*fk3_4;
u_5 = u_1+(439/216)*fk4_1 - 8*fk4_2 + (3680/513)*fk4_3 - (845/4104)*fk4_4;
theta_5 = theta_1+(439/216)*fk5_1 - 8*fk5_2 + (3680/513)*fk5_3 - (845/4104)*fk5_4;
fk3_5 = deltat*u_5;
fk4_5 = deltat*(-gm*r_5*r_5^(-gam)/(a+r_5)^(3-gam)+angmom2/(r_5^3.0));
fk5_5 = deltat*(angmom/(r_5^2.0));
r_6 = r_1-(8/27)*fk3_1 - 2*fk3_2 - (3544/2565)*fk3_3 + (1859/4104)*fk3_4-(11/40)*fk3_5;
u_6 = u_1-(8/27)*fk4_1 - 2*fk4_2 - (3544/2565)*fk4_3 + (1859/4104)*fk4_4-(11/40)*fk4_5;
theta_6 = theta_1-(8/27)*fk5_1 - 2*fk5_2 - (3544/2565)*fk5_3 + (1859/4104)*fk5_4-(11/40)*fk5_5;
fk3_6 = deltat*u_6;
fk4_6 = deltat*(-gm*r_6*r_6^(-gam)/(a+r_6)^(3-gam)+angmom2/(r_6^3.0));
fk5_6 = deltat*(angmom/(r_6^2.0));
fm3_1 = m1 + 25*fk3_1/216+1408*fk3_3/2565+2197*fk3_4/4104-fk3_5/5;
fm4_1 = m2 + 25*fk4_1/216+1408*fk4_3/2565+2197*fk4_4/4104-fk4_5/5;
fm5_1 = m3 + 25*fk5_1/216+1408*fk5_3/2565+2197*fk5_4/4104-fk5_5/5;
fm3_2 = m1 + 16*fk3_1/135+6656*fk3_3/12825+28561*fk3_4/56430-9*fk3_5/50+2*fk3_6/55;
fm4_2 = m2 + 16*fk4_1/135+6656*fk4_3/12825+28561*fk4_4/56430-9*fk4_5/50+2*fk4_6/55;
fm5_2 = m3 + 16*fk5_1/135+6656*fk5_3/12825+28561*fk5_4/56430-9*fk5_5/50+2*fk5_6/55;
R3 = abs(fm3_1-fm3_2)/deltat;
R4 = abs(fm4_1-fm4_2)/deltat;
R5 = abs(fm5_1-fm5_2)/deltat;
err3 = 0.84*(epsilon/R3)^(1/4);
err4 = 0.84*(epsilon/R4)^(1/4);
err5 = 0.84*(epsilon/R5)^(1/4);
if R3<= epsilon
time = time+deltat;
fm3 = fm3_1;
i = i+1;
deltat = err3*deltat;
end
if R4<= epsilon
time = time+deltat;
fm4 = fm4_1;
i = i+1;
deltat = err4*deltat;
end
if R5<= epsilon
time = time+deltat;
fm5 = fm5_1;
i = i+1;
deltat = err5*deltat;
end
e=2*gm^2.0/(2*angmom2);
ecc=(1.0+(2.0*e*angmom2)/(gm^2.0))^0.5;
x(i)=r_1*cos(theta_1)/(1000.0*parsec);
y(i)=r_1*sin(theta_1)/(1000.0*parsec);
time=time+deltat;
r(i)=r_1;
time1(i)=time;
end
figure()
plot(x,y, '-k');
TeXString = title('Plot of Orbit in Gamma Potential Obtained Using RK4')
xlabel('x')
ylabel('y')
2 ответа
Вы получаете сложные значения, потому что в какой-то момент npts - time < 0, Вы можете распечатать значения deltat проверить ошибку.
Кроме того, ваш код, по-видимому, не учитывает случай, когда оценка ошибки превышает допустимое отклонение. Когда ваша оценка ошибки превышает допустимое отклонение, вы должны:
- Перенесите время и решение
- рассчитать новый размер шага на основе формулы, и
- Пересчитайте ваше решение и оценку ошибки.
Тот факт, что вы не знаете, сколько итераций вам придется пройти, делает цикл for для адаптивной ступени Kutta немного неловким. Вместо этого я предлагаю использовать цикл while.
Вы используете "я" в своем коде. "i" возвращает базовую мнимую единицу. "i" эквивалентно sqrt(-1). Попробуйте использовать другой идентификатор в ваших циклах и используйте только "i" или "j" в вычислениях, где используются комплексные числа.