用ODE45在matlab中求解一个四阶ODE系统

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

这不是一个完整的解决方案，但比评论更长。

您似乎将t theta和phi分别用于两个目的。例如，可能会更改时间刻度行以使用s。

看看你的rhs，我认为它应该有两个参数。第一个参数应该是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