Решение системы из 4 ODE второго порядка в Matlab с использованием ODE45

Question

Мне нужно решить эту систему уравнений второго порядка, используя ODE45 в Matlab. Я знаком только с использованием ODE45, возможно, с одним или двумя уравнениями, но не так много. Вот что у меня есть, но я не уверен, как это исправить:

function second_oder_ode

t = 0:0.001:3;   % time scale
theta = pi/2;
phi = 0;

initial_t    = 0;
initial_dtds = 0;
initial_r    = 0;
initial_drds = 0;
initial_theta    = 0;
initial_dthetads = 0;
initial_phi    = 0;
initial_dphids = 0;

[t,x] = ode45( @(t,y) rhs(t,y,theta,phi), t, [initial_t initial_dtds initial_r initial_drds initial_theta initial_dthetads initial_phi initial_dphids] );

plot(t,x(:,1));
xlabel('t'); ylabel('r');
%%%%%%%%%%%%%% STATE VECTORS %%%%%%%%%%%%%%%%%%%%%
t = [t, dtds];
r = [r, drds];
theta = [theta, dthetads];
phi = [phi, dphids];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dxds=rhs(t,r, theta, phi)
    dxds_1 = t(2);
    dxds_2 = -((2/3)*R0^(2)*((3*H0/(2*cv))^(4/3))*y(1)^(1/3))*(y(4)^(2)+abs(y(3))^(2)*y(6)^(2)+abs(y(3))^(2)*(sin(y(5)))^(2)*y(8)^(2));
    dxds_3 = r(2);
    dxds_4 = -(4/(3*y(1)))*y(2)*y(4)+abs(y(3))*(y(6)^(2)+(sin(y(5)))^(2)*y(8)^(2));
    dxds_5 = theta(2);
    dxds_6 = -(2/abs(y(3)))*y(4)*y(6)-(4/(3*y(1)))*y(2)*y(6)+((sin(2*y(5)))*(y(8)^(2))/2);
    dxds_7 = phi(2);
    dxds_8 = -2*y(8)*((2/(3*y(1)))*y(2)+(cot(y(5)))*y(6)+(1/abs(y(3)))*y(4));

    dxdt=[dxdt_1; dxdt_2; dxdt_3; dxdt_4; dxdt_5; dxdt_6; dxdt_7; dxdt_8];

   end
end

Answer 1

Это не полное решение, но оно длиннее комментария.

Вы, кажется, используете ttheta а также phi каждый для двух целей. Например, возможно, измените шкалу времени, чтобы использовать s,

Глядя на ваш одефун, тут звонили rhsЯ думаю, это должно иметь 2 аргумента. Первый аргумент должен быть s (даже если он не используется) и второй аргумент должен быть векторным вводом, где значения y(1:8) соответствовать [t; t_s; r; r_s; theta; theta_s; phi; phi_s],

Выражения для dxds_# должен содержать только условия y и нет t, r, theta или же phi,

У вас есть опечатки на линии dxdt= который должен содержать dxds термины.

Ниже я сделал некоторые изменения. Я не проверял это; У меня нет констант R0, H0 или же cv и я думаю, что это будет равномерно ноль.

Также я удалил первые вхождения theta а также phi, Я не уверен, как они должны быть представлены в rhs,

function second_oder_ode

s = 0:0.001:3;   % time scale
%theta = pi/2;
%phi = 0;

initial_t    = 0;
initial_dtds = 0;
initial_r    = 0;
initial_drds = 0;
initial_theta    = 0;
initial_dthetads = 0;
initial_phi    = 0;
initial_dphids = 0;

y0 = [initial_t, initial_dtds, initial_r, initial_drds, ...
      initial_theta, initial_dthetads, initial_phi, initial_dphids];

[s,x] = ode45( @(s,y) rhs(s,y), s, y0);

plot(s,x(:,1));
xlabel('s'); ylabel('t(s)');
%%%%%%%%%%%%%% STATE VECTORS %%%%%%%%%%%%%%%%%%%%%
%t = [t, dtds];
%r = [r, drds];
%theta = [theta, dthetads];
%phi = [phi, dphids];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dxds=rhs(s,y)
    dxds_1 = y(2);
    dxds_2 = -((2/3)*R0^(2)*((3*H0/(2*cv))^(4/3))*y(1)^(1/3))*(y(4)^(2)+abs(y(3))^(2)*y(6)^(2)+abs(y(3))^(2)*(sin(y(5)))^(2)*y(8)^(2));
    dxds_3 = y(4);
    dxds_4 = -(4/(3*y(1)))*y(2)*y(4)+abs(y(3))*(y(6)^(2)+(sin(y(5)))^(2)*y(8)^(2));
    dxds_5 = y(6);
    dxds_6 = -(2/abs(y(3)))*y(4)*y(6)-(4/(3*y(1)))*y(2)*y(6)+((sin(2*y(5)))*(y(8)^(2))/2);
    dxds_7 = y(8);
    dxds_8 = -2*y(8)*((2/(3*y(1)))*y(2)+(cot(y(5)))*y(6)+(1/abs(y(3)))*y(4));

    dxds=[dxds_1; dxds_2; dxds_3; dxds_4; dxds_5; dxds_6; dxds_7; dxds_8];

   end
end