Reducing run time for vpasolve

Hello, I am a very poor coder and need to create a program that solves for the concentration of 3 different chemical species in a set of CSTR reactors either in series or parallel. The concentration depends on the residence time of the reactor. To do this I am using vpasolve to solve my set of kinetic equations, however my program takes too long to run to completetion (not sure on the exact time but I have let it run for 30+ minutes and it has yet to produce a graph. How can I reduce the time it takes to solve the set of equations? Attached is my MATLAB file and below is the text of my code:

function main
CSTRseries
CSTRparallel
end


function CSTRseries

%variables
global k1 k2 k3 k4 Cao Cro Cso
k1=.002777778; %Rate constant (mol/(L*s))
k2=.0002777778; %Rate constant (mol/(L*s))
k3=.0002777778; %Rate constant (mol/(L*s))
k4=.0002777778; %Rate constant (mol/(L*s))
Cao=5; %Initial concentration of A (mol/L)
Cro=0; %Initial concentration of R (mol/L)
Cso=0; %Initial concentration of S (mol/L)

%Solving using the numerical ode solver (numerical solution)
syms Ca Cr Cs
tau=zeros(1,3600);
C=zeros(3,3600);
for i=1:3600
    tau(i)=i;
    [x,y,z] = vpasolve([(Ca-Cao)==(-k1*Ca + k2*Cr)*tau(i),(Cr-Cro)==(k1*Ca + -k2*Cr + -k3*Cr + k4*Cs)*tau(i),(Cs-Cso)==(k3*Cr + -k4*Cs)*tau(i)],[Ca,Cr,Cs],[0.0000001 Inf;0.00000001 Inf;0.0000001 Inf]);
    C(1,i)=x;
    C(2,i)=y;
    C(3,i)=z;
end

%Plot data
figure(1)
plot(tau,C)
title('CSTRs in Series')
xlabel('tau')
ylabel('Concentration')
legend('A','R','S')

end

function CSTRparallel

%variables
global k1 k2 k3 k4 Cao Cro Cso
k1=.000833333; %Rate constant (mol/(L*s))
k2=.000833333; %Rate constant (mol/(L*s))
k3=.0002777778; %Rate constant (mol/(L*s))
k4=.0002777778; %Rate constant (mol/(L*s))
Cao=5; %Initial concentration of A (mol/L)
Cro=0; %Initial concentration of R (mol/L)
Cso=0; %Initial concentration of S (mol/L)

%Solving using the numerical ode solver (numerical solution)
syms Ca Cr Cs
tau=zeros(1,3600);
C=zeros(3,3600);
for i=1:3600
    [x,y,z] = vpasolve([(Ca-Cao)==(-k1*Ca + k2*Cr + -k3*Ca + k4*Cs)*tau(i),(Cr-Cro)==(k1*Ca + -k2*Cr)*tau(i),(Cs-Cso)==(k3*Ca + -k4*Cs)*tau(i)],[Ca,Cr,Cs],[0.0000001 5;0.00000001 5;0.0000001 5]);
    C(1,i)=x;
    C(2,i)=y;
    C(3,i)=z;
end

%Plot data
figure(2)
plot(tau,C)
title('CSTRs in Parallel')
xlabel('tau')
ylabel('Concentration')
legend('A','R','S')

end

[Ameer Hamza]

Is there a specific reason for not using a numeric solver, like fsolve?

[Samuel Dixon]

I believe our professor wants us to use vpasolve, as this was part of an example she provided.

[Ameer Hamza]

Are you solving exactly same equation 3600 times?

[x,y,z] = vpasolve([(Ca-Cao)==(-k1*Ca + k2*Cr + -k3*Ca + k4*Cs)*tau(i),(Cr-Cro)==(k1*Ca + -k2*Cr)*tau(i),(Cs-Cso)==(k3*Ca + -k4*Cs)*tau(i)],[Ca,Cr,Cs],[0.0000001 5;0.00000001 5;0.0000001 5]);

tau(i) are all zeros so nothing change in each iteration.

[Samuel Dixon]

I believe thats what my problem is, but I am unsure how to fix it. Tau should be [1,2,3...3600] not an array of all zeros. But I do not know how to prevent MATLAB from solving the entire array again every time it moves to the next iterations. Basically I want it to solve for Ca at each value of tau from 1 to 3600 and then graph those Ca values against the corresponding tau value.

Answer 1

[Ameer Hamza]

You can get some speed improvement by using the following method. I pre-solved the equation in terms of tau variable, and inside the for loop, I substituted the value of tau. I changed one of the functions and ran for tau_=1:100 (100 elements). I observed a speed improvement of about 2x. You can similarly change other function following this method.

function CSTRseries

%variables
global k1 k2 k3 k4 Cao Cro Cso
k1=.002777778; %Rate constant (mol/(L*s))
k2=.0002777778; %Rate constant (mol/(L*s))
k3=.0002777778; %Rate constant (mol/(L*s))
k4=.0002777778; %Rate constant (mol/(L*s))
Cao=5; %Initial concentration of A (mol/L)
Cro=0; %Initial concentration of R (mol/L)
Cso=0; %Initial concentration of S (mol/L)

%Solving using the numerical ode solver (numerical solution)
syms Ca Cr Cs tau
sol = solve([(Ca-Cao)==(-k1*Ca + k2*Cr)*tau, (Cr-Cro)==(k1*Ca + -k2*Cr + -k3*Cr + k4*Cs)*tau,(Cs-Cso)==(k3*Cr + -k4*Cs)*tau], [Ca,Cr,Cs]);
tau_=1:100;
C=zeros(3,100);
for i=1:100
    C(1,i) = subs(sol.Ca, tau, tau_(i));
    C(2,i) = subs(sol.Cr, tau, tau_(i));
    C(3,i) = subs(sol.Cs, tau, tau_(i));
end

%Plot data
figure(1)
plot(tau_,C)
title('CSTRs in Series')
xlabel('tau')
ylabel('Concentration')
legend('A','R','S')

end

[Samuel Dixon]

Im getting the error:

"Error using syms (line 255)

Unable to create a symbolic variable 'tau' by using 'syms' inside a MATLAB function because 'tau' is an existing function name, class name, method name, and so

on. Use 'sym' instead.

Error in CSTRtest (line 12)

syms Ca Cr Cs tau"

Is this because tau is both a symbol and array?

[Ameer Hamza]

Yes, in my code, I changed the numeric array name from tau to tau_

[Samuel Dixon]

I copied your code exactly from your comment and I am still getting the same error.

[Ameer Hamza]

Can you paste you new code here? I guess there was some mistake in naming the variables.

[Samuel Dixon]

I managed to figure it out, although I'm not quite sure what I did to fix it either. Thank you for your help though!

[Ameer Hamza]

Glad to be of help.