MATLAB Answers

Simultaneous differential equations - derivative of y^3 wrt t

2 views (last 30 days)
Nora Rafael
Nora Rafael on 16 Nov 2019
Answered: darova on 29 Nov 2019
I have been trying to solve this problem for sometime and so far I have received an answer that requires the symbolic toolbox which I dont have - I would appreciate any help!
I need to solve these 2 differential equations simultaneously.
cube derivative.JPG
But I dont know how to code the dr*^3/dt. My code is below:
function dydt=odefcnNY_v3(t,y,D,Cs,rho,r0,N,V,Af)
dydt=zeros(2,1);
dydt(1)=(-3*D*Cs/rho*r0^2)*y(1)*(1-y(2));
dydt(2)=(D*4*pi*N*r0*(1-y(2))*y(1)-(Af*y(2)))/V;
end
y(1) = r* and
y(2) = C*
So
dydt(1) = dr*/dt and
dydt(2) = dC*/dt
In my case dydt(1) needs to be replaced with something that would solve dr*^3/dt and not dr*/dt. The rest of the code is below.
D=4e-9;%m2/s
rho=1300; %kg/m3
r0=10.1e-6; %m dv50
Cs=0.0016; %kg/m3
V=1.5e-6;%m3
W=4.5e-8; %kg
N=W/(4/3*pi*r0^3*rho);
Af=0.7e-6/60; %m3/s
tspan=[0 24*3600]; %s in 24 hrs
y0=[r0 0];
[t,y]=ode45(@(t,y) odefcnNY_v3(t,y,D,Cs,rho,r0,Af,N,V), tspan, y0);
plot(t/3600,y(:,1),'-o') %plot time in hr, and r*
xlabel('time, hr')
ylabel('radius,um')
legend('DCU')
plot(t/3600,y(:,2),'-') %plot time in hr, and C*
xlabel('time,hr')
ylabel('C* (C/Cs)')
legend('DCU')

  9 Comments

Show 6 older comments
David Goodmanson
David Goodmanson on 21 Nov 2019
There appears to be a mistake on the second slide. For dC*/dt the highlighted equation has a term Af*(C_L), whereas your equation uses Af*(C*). But if you look at the derivation on that slide, Af*(C*) would be correct.
Nora Rafael
Nora Rafael on 29 Nov 2019
Hi David, that's right, the slide is wrong - it should be C*. I originally had the equations in terms of CL and later converted to C*, and didn't update this part of the equation on the slides.
Nora Rafael
Nora Rafael on 29 Nov 2019
I am now trying to implement a similar problem, but with slightly different parameters and a simplified function with no AfxC* term. This one seems to be taking hours and still no solution, whereas the previous one was taking only a few seconds.

Sign in to comment.

Answers (1)

darova
darova on 29 Nov 2019
Just get
dydt(1)=(-3*D*Cs/rho/r0^2)*y(1)*(1-y(2));
dydt(1)=nthroots(dydt(1),3);

  0 Comments

Sign in to comment.

Sign in to answer this question.