Derivative of y^3 using Matlab ode45

Question

I have two differential equations I would like to solve simultaneously using ode45 in Matlab, shown below. The problem I have is that the second differential equation is actually d(y^3)/dt, i.e. the derivative of y^3 with respect to t. How do I express this?

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),'-.')

jman jman · Accepted Answer · 2019-11-01T22:33:31

d(y^3)/dt = 3 * y^2 * dy/dt so you can divide the second equation with (3 * y(2) * y(2)) . This should be fine unless y(2) gets close to zero.

If it does get close to zero then you will have to re-write your differential equation in terms of another variable in which the second equation doesn't have singularities. For example you can have u(t) = y^3(t) and then y(t) = u(t)^{1/3} (you will have to figure out if there will be ambiguities with negative cube roots).

Derivative of y^3 using Matlab ode45

1 Answers