How can I solve these differential equations?

Hi,

I am trying to solve these three differential equations using the explicit equations I have. I also want to plot Fa vs. V, Fb vs. V, and Fc vs. V as shown in Figure E6-11. Please help me.

Thank you.

[James Tursa]

What have you tried so far? I assume you want a numerical solution using ode45(), so that would mean a 3-element state vector representing [Fa,Fb,Fc] with independent variable V. Have you looked into the ode45() examples, particularly the ones with multiple variables?

I would note it is suspicious that the following variables appear to be unused in the equations: Fao, vo, Tau, X, rateA

[Huan Gu]

Hi James,

Thank you for responding. I have tried to make dFa/dv=a(1-Fa) work, but if I try to use the same code to solve dFa/dv=a*((1-Fa)/Fa)^2, I had a lot of errors that I cannot fix.

H.

[Huan Gu]

If I combine the explicit information, d(Fa)/dV=-k*(Ct0^2)*((Fa)^2)/((3/2*(Fa0)-1/2*(Fa))^2)

here: k=274.4

Ft=Fa+Fb+Fc

Fa=Fa0-Fa0*X

Fb=Fa0*X

Fc=1/2*Fa0*X

Ct0=1641/8.314/T

T=698

Fa0=0.00000226

[Sam Chak]

James was right. Items #7, 8, 9, 11 are unused. I listed out , , , and then performed the back substitutions.

E = 24000;
T = 697;
k = 0.29*exp(E/1.987*(1/500 - 1/T));
Cto = 1641/8.314/T;
Ca = Cto*y(1)/(y(1) + y(2) + y(3));
ra = - k*Ca^2;
rb = -ra;
rc = -ra/2;

Answer 1

[James Tursa]

A general outline would be to create your derivative function as follows:

function dy = myderivative(t,y,T,E)
Fa = y(1);
Fb = y(2);
Fc = y(3);
ra = _____; % you fill this in
rb = _____; % you fill this in
rc = _____; % you fill this in
dy = [ra;rb;rc];
end

Then in your main code have a function handle for this:

T = 698;
E = 24000;
f = @(t,y) myderivative(t,y,T,E)

Then you can use f along with some initial conditions and V span to call ode45().

Answer 2

[Sam Chak]

Hi Huan Gu

I've got these results, not exactly the same patterns, but very close.

clear all; clc
Vspan = 0:1e-9:1e-5;
y0 = [2.26e-5; 0; 0];
[V, y] = ode45(@(t,y) [- (0.29*exp(24000/1.987*(1/500 - 1/697)))*((1641/8.314/697)*y(1)/(y(1) + y(2) + y(3)))^2; (0.29*exp(24000/1.987*(1/500 - 1/697)))*((1641/8.314/697)*y(1)/(y(1) + y(2) + y(3)))^2; 0.5*(0.29*exp(24000/1.987*(1/500 - 1/697)))*((1641/8.314/697)*y(1)/(y(1) + y(2) + y(3)))^2], 
Vspan, y0);
plot(V, y(:,1), 'linewidth', 2, 'Color', [0.4940, 0.1840, 0.5560], V, y(:,2), 'linewidth', 2, 'Color', [0.4660, 0.6740, 0.1880], V, y(:,3), 'linewidth', 2, 'Color', [0.9290, 0.6940, 0.1250])
legend('Fa', 'Fb', 'Fc')
grid on
xlabel('V (dm^3)')
ylabel('Flow rates (mol/s)')