MATLAB: Can’t figure out how to finalize the code

matlab functionMATLAB Online Server

I need to write a code in order to get Ca1 Ca2 Ca3 Ca4 Cb1 Cb2 Cb3 Cb4 by using this part which is given below. Only for "Eqns_2nd_Order". But I didn't figure out how can I do. So I need to help, could you help me please?

Eqns_2nd_Order=[(Qin*Cain)-(Q12*Ca1)-(k21*(Cain)^2*V1)%=0























(Q12*Ca1)+(Q32*Ca3)-(Q23*Ca2)-(k22*(Ca2)^2*(V2))%=0
(Q23*Ca2)+(Q43*Ca4)-(Q32*Ca3)-(Q34*Ca3)-(k23*(Ca3)^2*(V3))%=0
(Q34*Ca3)-(Q43*Ca4)-(Qout*Ca4)-(k24*(Ca4)^2*V4)%=0
(Qin*Cbin)-(Q12*Cb1)+(k21*(Cain)^2*V1)%=0
(Q12*Cb1)+(Q32*Cb3)-(Q23*Cb2)+(k22*(Ca2)^2*(V2))%=0
(Q23*Cb2)+(Q43*Cb4)-(Q32*Cb3)-(Q34*Cb3)+(k23*(Ca3)^2*(V3))%=0
(Q34*Cb3)-(Q43*Cb4)-(Qout*Cb4)+(k24*(Ca4)^2*V4)]%=0
%Cain=initial concentration of chemical A in units of ML^-3
%Cbin=initial concentration of chemical B in units of ML^-3

MAIN CODE:


%V1=volume of first reactor in L
%V2=volume of second reactor in L
%V3=volume of third reactor in L
%V4=volume of forth reactor in L
 
%k1=first order reaction rate in units of (L^3)(M^-1)(T^-1)
%k2=second order reaction rate in units of (L^3)(M^-1)(T^-1)
%Each reactor has different reaction rate
 
%M=mass in units of M
%Qin=inflow in units of (L^3)(T^-1)
%Qout=outflow in units of (L^3)(T^-1)
%Qxy=the flow from reactor x to reactor y (x&y can be 1,2,3..)
%UNITS
%The unit for V1, V2, V3 & V4 is L
%The unit for k1 is hr^-1
%The unit for k2 is L/(mg*hr)
%The unit for M is mg
%The unit for Qin, Qout,Q12, Q23, Q34, Q32, Q43 & Qxy is L/hr
%Given and calculated values  
Qin=10;
Q12=10;
Q23=10+5; %Q23=Q12+Q32
Q34=15-5+3; %Q34=Q23-Q32+Q43
Qout=10;
Q32=5;
Q43=3;
 
Cain=1;
Cbin=0;
 
V1=25;
V2=75;
V3=100;
V4=25;
 
k11=0.05;
k12=0.1;
k13=0.5;
k14=0.1;
 
k21=0.04;
k22=0.08;
k23=0.4;
k24=0.08;
 
%% PARTA
 
%Mass Balance Eqns for Chemical Species A
 
%For reactor 1

%(Qin*Cain)-(Q12*Ca1)-(k11*Cain*V1)==0  %for first order reaction
%(Qin*Cain)-(Q12*Ca1)-(k21*(Cain)^2*V1)==0  %for second oreder reaction
 
%For reactor 2

%(Q12*Ca1)+(Q32*Ca3)-(Q23*Ca2)-(k12*Ca2*V2)==0 %for first order reaction
%(Q12*Ca1)+(Q32*Ca3)-(Q23*Ca2)-(k22*(Ca2)^2*(V2))==0 %for second oreder reaction
 
%For reactor 3

%(Q23*Ca2)+(Q43*Ca4)-(Q32*Ca3)-(Q34*Ca3)-(k13*Ca3*V3)==0  %for first order reaction
%(Q23*Ca2)+(Q43*Ca4)-(Q32*Ca3)-(Q34*Ca3)-(k23*(Ca3)^2*(V3))==0  %for second order reaction
 
%For reactor 4

%(Q34*Ca3)-(Q43*Ca4)-(Qout*Ca4)-(k14*Ca4*V4)==0  %for first order reaction
%(Q34*Ca,3)-(Q43*Ca,4)-(Qout*Ca,4)-(k24*(Ca4)^2*V4)==0  %for second order reaction
 
%Mass Balance Eqns for Chemical Species B
 
%For reactor 1
%(Qin*Cbin)-(Q12*Cb1)+(k11*Ca1*V1)==0  %for first order reaction
%(Qin*Cbin)-(Q12*Cb1)+(k21*(Ca1)^2*V1)==0  %for second oreder reaction
 
%For reactor 2
%(Q12*Cb1)+(Q32*Cb3)-(Q23*Cb2)+(k12*Ca2*V2)==0 %for first order reaction
%(Q12*Cb1)+(Q32*Cb3)-(Q23*Cb2)+(k22*(Ca2)^2*(V2))==0 %for second oreder reaction
 
%For reactor 3
%(Q23*Cb2)+(Q43*Cb4)-(Q32*Cb3)-(Q34*Cb3)+(k13*Ca3*V3)==0  %for first order reaction
%(Q23*Cb2)+(Q43*Cb4)-(Q32*Cb3)-(Q34*Cb3)+(k23*(Ca3)^2*(V3))==0  %for second order reaction
 
%For reactor 4
%(Q34*Cb3)-(Q43*Cb4)-(Qout*Cb4)+(k14*Ca4*V4)==0  %for first order reaction
%(Q34*Cb3)-(Q43*Cb4)-(Qout*Cb4)+(k24*(Ca4)^2*V4)==0  %for second order reaction
 
syms Ca1 Ca2 Ca3 Ca4 Cb1 Cb1 Cb2 Cb3 Cb4
Eqns_1st_Order=[(Qin*Cain)-(Q12*Ca1)-(k11*Cain*V1)%=0
(Q12*Ca1)+(Q32*Ca3)-(Q23*Ca2)-(k12*Ca2*V2)%=0
(Q23*Ca2)+(Q43*Ca4)-(Q32*Ca3)-(Q34*Ca3)-(k13*Ca3*V3)%=0
(Q34*Ca3)-(Q43*Ca4)-(Qout*Ca4)-(k14*Ca4*V4)%=0
(Qin*Cbin)-(Q12*Cb1)+(k11*Cain*V1)%=0
(Q12*Cb1)+(Q32*Cb3)-(Q23*Cb2)+(k12*Ca2*V2)%=0
(Q23*Cb2)+(Q43*Cb4)-(Q32*Cb3)-(Q34*Cb3)+(k13*Ca3*V3)%=0
(Q34*Cb3)-(Q43*Cb4)-(Qout*Cb4)+(k14*Ca4*V4)]%=0
 
Eqns_2nd_Order=[(Qin*Cain)-(Q12*Ca1)-(k21*(Cain)^2*V1)%=0
(Q12*Ca1)+(Q32*Ca3)-(Q23*Ca2)-(k22*(Ca2)^2*(V2))%=0
(Q23*Ca2)+(Q43*Ca4)-(Q32*Ca3)-(Q34*Ca3)-(k23*(Ca3)^2*(V3))%=0
(Q34*Ca3)-(Q43*Ca4)-(Qout*Ca4)-(k24*(Ca4)^2*V4)%=0
(Qin*Cbin)-(Q12*Cb1)+(k21*(Cain)^2*V1)%=0
(Q12*Cb1)+(Q32*Cb3)-(Q23*Cb2)+(k22*(Ca2)^2*(V2))%=0
(Q23*Cb2)+(Q43*Cb4)-(Q32*Cb3)-(Q34*Cb3)+(k23*(Ca3)^2*(V3))%=0
(Q34*Cb3)-(Q43*Cb4)-(Qout*Cb4)+(k24*(Ca4)^2*V4)]%=0

Best Answer

%V1=volume of first reactor in L
%V2=volume of second reactor in L
%V3=volume of third reactor in L
%V4=volume of forth reactor in L
 
%k1=first order reaction rate in units of (L^3)(M^-1)(T^-1)
%k2=second order reaction rate in units of (L^3)(M^-1)(T^-1)
%Each reactor has different reaction rate
 
%M=mass in units of M
%Qin=inflow in units of (L^3)(T^-1)
%Qout=outflow in units of (L^3)(T^-1)
%Qxy=the flow from reactor x to reactor y (x&y can be 1,2,3..)
%UNITS
%The unit for V1, V2, V3 & V4 is L
%The unit for k1 is hr^-1
%The unit for k2 is L/(mg*hr)
%The unit for M is mg
%The unit for Qin, Qout,Q12, Q23, Q34, Q32, Q43 & Qxy is L/hr
%Given and calculated values  
Qin=10;
Q12=10;
Q23=10+5; %Q23=Q12+Q32
Q34=15-5+3; %Q34=Q23-Q32+Q43
Qout=10;
Q32=5;
Q43=3;
 
Cain=1;
Cbin=0;
 
V1=25;
V2=75;
V3=100;
V4=25;
 
k11=0.05;
k12=0.1;
k13=0.5;
k14=0.1;
 
k21=0.04;
k22=0.08;
k23=0.4;
k24=0.08;
 
%% PARTA
 
%Mass Balance Eqns for Chemical Species A
 
%For reactor 1

%(Qin*Cain)-(Q12*Ca1)-(k11*Cain*V1)==0  %for first order reaction
%(Qin*Cain)-(Q12*Ca1)-(k21*(Cain)^2*V1)==0  %for second oreder reaction
 
%For reactor 2

%(Q12*Ca1)+(Q32*Ca3)-(Q23*Ca2)-(k12*Ca2*V2)==0 %for first order reaction
%(Q12*Ca1)+(Q32*Ca3)-(Q23*Ca2)-(k22*(Ca2)^2*(V2))==0 %for second oreder reaction
 
%For reactor 3

%(Q23*Ca2)+(Q43*Ca4)-(Q32*Ca3)-(Q34*Ca3)-(k13*Ca3*V3)==0  %for first order reaction
%(Q23*Ca2)+(Q43*Ca4)-(Q32*Ca3)-(Q34*Ca3)-(k23*(Ca3)^2*(V3))==0  %for second order reaction
 
%For reactor 4

%(Q34*Ca3)-(Q43*Ca4)-(Qout*Ca4)-(k14*Ca4*V4)==0  %for first order reaction
%(Q34*Ca,3)-(Q43*Ca,4)-(Qout*Ca,4)-(k24*(Ca4)^2*V4)==0  %for second order reaction
 
%Mass Balance Eqns for Chemical Species B
 
%For reactor 1
%(Qin*Cbin)-(Q12*Cb1)+(k11*Ca1*V1)==0  %for first order reaction
%(Qin*Cbin)-(Q12*Cb1)+(k21*(Ca1)^2*V1)==0  %for second oreder reaction
 
%For reactor 2
%(Q12*Cb1)+(Q32*Cb3)-(Q23*Cb2)+(k12*Ca2*V2)==0 %for first order reaction
%(Q12*Cb1)+(Q32*Cb3)-(Q23*Cb2)+(k22*(Ca2)^2*(V2))==0 %for second oreder reaction
 
%For reactor 3
%(Q23*Cb2)+(Q43*Cb4)-(Q32*Cb3)-(Q34*Cb3)+(k13*Ca3*V3)==0  %for first order reaction
%(Q23*Cb2)+(Q43*Cb4)-(Q32*Cb3)-(Q34*Cb3)+(k23*(Ca3)^2*(V3))==0  %for second order reaction
 
%For reactor 4
%(Q34*Cb3)-(Q43*Cb4)-(Qout*Cb4)+(k14*Ca4*V4)==0  %for first order reaction
%(Q34*Cb3)-(Q43*Cb4)-(Qout*Cb4)+(k24*(Ca4)^2*V4)==0  %for second order reaction
 
syms Ca1 Ca2 Ca3 Ca4 Cb1 Cb2 Cb3 Cb4
Eqns_1st_Order=[(Qin*Cain)-(Q12*Ca1)-(k11*Cain*V1)%=0















(Q12*Ca1)+(Q32*Ca3)-(Q23*Ca2)-(k12*Ca2*V2)%=0
(Q23*Ca2)+(Q43*Ca4)-(Q32*Ca3)-(Q34*Ca3)-(k13*Ca3*V3)%=0
(Q34*Ca3)-(Q43*Ca4)-(Qout*Ca4)-(k14*Ca4*V4)%=0
(Qin*Cbin)-(Q12*Cb1)+(k11*Cain*V1)%=0
(Q12*Cb1)+(Q32*Cb3)-(Q23*Cb2)+(k12*Ca2*V2)%=0
(Q23*Cb2)+(Q43*Cb4)-(Q32*Cb3)-(Q34*Cb3)+(k13*Ca3*V3)%=0
(Q34*Cb3)-(Q43*Cb4)-(Qout*Cb4)+(k14*Ca4*V4)]%=0
Eqns_1st_Order = 
 
Eqns_2nd_Order=[(Qin*Cain)-(Q12*Ca1)-(k21*(Cain)^2*V1)%=0
(Q12*Ca1)+(Q32*Ca3)-(Q23*Ca2)-(k22*(Ca2)^2*(V2))%=0
(Q23*Ca2)+(Q43*Ca4)-(Q32*Ca3)-(Q34*Ca3)-(k23*(Ca3)^2*(V3))%=0
(Q34*Ca3)-(Q43*Ca4)-(Qout*Ca4)-(k24*(Ca4)^2*V4)%=0
(Qin*Cbin)-(Q12*Cb1)+(k21*(Cain)^2*V1)%=0
(Q12*Cb1)+(Q32*Cb3)-(Q23*Cb2)+(k22*(Ca2)^2*(V2))%=0
(Q23*Cb2)+(Q43*Cb4)-(Q32*Cb3)-(Q34*Cb3)+(k23*(Ca3)^2*(V3))%=0
(Q34*Cb3)-(Q43*Cb4)-(Qout*Cb4)+(k24*(Ca4)^2*V4)]%=0
Eqns_2nd_Order = 
sol = solve(Eqns_2nd_Order, [Ca1 Ca2 Ca3 Ca4 Cb1 Cb2 Cb3 Cb4])
sol = struct with fields:
    Ca1: [8×1 sym]
    Ca2: [8×1 sym]
    Ca3: [8×1 sym]
    Ca4: [8×1 sym]
    Cb1: [8×1 sym]
    Cb2: [8×1 sym]
    Cb3: [8×1 sym]
    Cb4: [8×1 sym]
sol.Ca1
ans = 
sol.Ca2
ans = 
vpa(sol.Ca2, 3)
ans =

Related Solutions

MATLAB: How to run script iteratively

 clc; clear; close all;
          format long 
          Ma=.84
          Ra=.287
          Rg=.287
          siga=1.4
          sigg=1.3333
          Pa=54.05
          Ta=225.7
          OPR=2:3:200
          tsfc=cell(1,67) ;      %PRE-ALLOCATION
          for i = 1:numel(OPR)
          Va=Ma*sqrt(siga*Ra*Ta)
          Cpoa= (siga*Ra)./(siga - 1)
          TINT = 1100:100:1800
          PRC = sqrt(OPR(i))
          % Diffuser
          efd=.93
          To1 = Ta*(1+((siga-1)/2)*(Ma^2))
          Toa=To1
          Po1 = Pa*(1 + efd*((siga - 1)/2)*(Ma^2))^(siga/siga-1)
          Poa=Pa*(Toa/Ta)^(siga/siga-1)
          %Compressors  
          efc=.87
          %Low Pressure
          Po2=Po1.*(PRC)
          To2 = To1 + ((To1*((PRC.^((siga-1)/siga)-1)))/efc)
          WC_LP = Cpoa*(To1-To2)
          %High Pressure
          To3 = To2 + ((To2.*((PRC.^((siga-1)/siga)-1)))/efc)
          Po3=Po2.*(PRC)
          WC_HP = Cpoa*(To2-To3)
          %Combustion
          efb =.98
          deltPob = .4 
          Po4= (1-deltPob)*Po3
          To4=TINT
          %Combustion Calculations for f
          MC=14.4
          MH=24.9
          MO=0
          Ycc=MC + (MH/4)- (MO/2)
          Tp=To4
          Tr=To3
          To= 298
          HRP_to = -8561991.6
          if Tp <= 1600
          hco2p = 56835 + 66.27*Tp + -11634.0*log(Tp)
          hh2op = 88923 + 49.36*Tp + -7940.8*log(Tp)
          hn2p= 31317 + 37.46*Tp + -4559.3*log(Tp)
          ho2p = 43388 + 42.27*Tp + -6635.4*log(Tp)
          else
          hco2p = 93048 + 68.58*Tp + -16979.0*log(Tp)
          hh2op = 154670 + 60.43*Tp + -19212.0*log(Tp)
          hn2p = 44639 + 39.32*Tp + -6753.4*log(Tp)
          ho2p =127010 + 46.25*Tp + -18798*log(Tp)   
          end
          if Tr <= 1600
          hco2r = 56835 + 66.27*Tr + -11634.0*log(Tr)
          hh2or = 88923 + 49.36*Tr + -7940.8*log(Tr)
          hn2r= 31317 + 37.46*Tr + -4559.3*log(Tr)
          ho2r = 43388 + 42.27*Tr + -6635.4*log(Tr)
          else
          hco2r = 93048 + 68.58*Tr + -16979.0*log(Tr)
          hh2or = 154670 + 60.43*Tr + -19212.0*log(Tr)
          hn2r = 44639 + 39.32*Tr + -6753.4*log(Tr)
          ho2r =127010 + 46.25*Tr + -18798*log(Tr)
          end
          hco2 = hco2p-hco2r
          hh2o = hh2op-hh2or
          hn2 = hn2p-hn2r
          ho2 =ho2p-ho2r
          h2oa=30.5
          h2ob=.0103
          co2a=28.8
          co2b=.028
          o2a=27
          o2b=.0079
          fuela= 38.1
          fuelb=.656
          Delt_HRP = (MH/2)*(h2oa*Tr +.5*h2ob*Tr.^2) - (h2oa*To +.5*h2ob*To.^2)+MC*(co2a*Tr +.5*co2b*Tr.^2) - (co2a*To +.5*co2b*To.^2)- Ycc*(o2a*Tr +.5*o2b*Tr.^2) - (o2a*To +.5*o2b*To.^2)-(fuela*Tr +.5*fuelb*Tr.^2) - (fuela*To +.5*fuelb*To.^2)
          HRP =HRP_to + Delt_HRP
          y = (-HRP - MC*hco2 - (MH/2).*hh2o + Ycc.*ho2)./(3.76.*hn2+ho2)
          fideal=197.7./(4.76*y*28.97)
          f=fideal*efb
          %Turbines | HIgh and Low have equal efficiencies 
          eft=.90
          efm=.99
          Cpog = (sigg*Rg)/(sigg-1)
          % High Pressure
          To5= To4 + (WC_HP./(efm*(1+f)*Cpog))
          Po5 = Po4*(1-((1-(To5/To4))/eft))^(sigg/(sigg-1))
          % Low Pressure
          To6= To5 + (WC_LP./(efm*(1+f)*Cpog))
          Po6 = Po5*(1-((1-(To6/To5))/eft))^(sigg/(sigg-1))
          %Nozzle
          %Choke Test
          efn = .95
          P_Po6 = (1-((1/efn)*(1-(2/(sigg+1)))))^(sigg/(sigg-1))
          R=Pa./Po6
          Me=1
          Pe=Po6*(P_Po6)
          Te=To6*(2/(sigg+1))
          roe= Pe./(Rg*Te)
          Ve= Me*sqrt(sigg*Rg*1000*Te)
          % Fs and T.S.F.C Calculations
          Fs = ((1+f).*Ve-Va) + (Pe-Pa).*((1+f)./(roe.*Ve))
          tsfc{i}= (f./Fs).*3600;
          end
          celldisp(tsfc)        %each subcell represents each value of OPR
tsfc = vertcat(tsfc{:})    %double array

MATLAB: Output of previous command into input

If those variables are in the workspace, it should run without any problems.

Ca1 = ...;
Ca2 = ...;
Ca3 = ...;
Ca4 = ...;
Cb1 = ...;
Cb2 = ...;
Cb3 = ...;
Cb4 = ...;
x0 = [Ca1 Ca2 Ca3 Ca4 Cb1 Cb2 Cb3 Cb4];
x = fsolve(@functiondescription, x0)

Best Answer

Related Solutions

MATLAB: How to run script iteratively

MATLAB: Output of previous command into input

Related Question