使用 ODE45 在 matlab 中求解 4 个二阶 ODE 的系统

Question

我需要在 matlab 中使用 ODE45 求解这个二阶方程组我只熟悉将 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.

查看您的 odefun，这里称为 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

使用 ODE45 在 matlab 中求解 4 个二阶 ODE 的系统

Solving a system of 4 second order ODE's in matlab using ODE45

matlab

ode45