matlab - 使用 Matlab ode45 对 y^3 求导

Question

我有两个微分方程，我想在 Matlab 中使用 ode45 同时求解，如下所示。我遇到的问题是二阶微分方程实际上是 d(y^3)/dt，即 y^3 对 t 的导数。我该如何表达？

function dydt=odefcnNY(t,y,D,Cs,rho,r0,Af,N,V)
dydt=zeros(2,1);
dydt(1)=(3*D*Cs/rho*r0^2)*(y(1)/r0)*(1-y(2)/Cs);
dydt(2)=(D*4*pi*r0*N*(1-y(2)/Cs)*(y(1)/r0)-(Af*y(2)/Cs))/V;
end

D=4e-9;
rho=1300;
r0=10.1e-6;
Cs=0.0016;
V=1.5e-6;
W=4.5e-8;
N=W/(4/3*pi*r0^3*rho);
Af=0.7e-6/60;
tspan=[0 75000];
y0=[r0 Cs];
[t,y]=ode45(@(t,y) odefcnNY(t,y,D,Cs,rho,r0,Af,N,V), tspan, y0);
plot(t,y(:,1),'-o',t,y(:,2),'-.')

Answer 1

d(y^3)/dt = 3 * y^2 * dy/dt所以你可以将第二个方程除以(3 * y(2) * y(2))。这应该没问题，除非 y(2) 接近于零。

如果它确实接近于零，那么您将不得不根据另一个变量重新编写您的微分方程，其中第二个方程没有奇点。例如你可以有u(t) = y^3(t)然后y(t) = u(t)^{1/3}（你必须弄清楚是否会有负立方根的歧义）。