%%%% solve ODE model to get output time profile % Program by Peter Halling (Uni Strathclyde), based on an original program by Hui Yu % Also calls ODE_ProchiralDiacid.m %clear all clear variables %following help guidance close all clc format long % Input values here may be changed showplot = false; S0 = 0.2; % M units B0 = 55; Q0 = 0; Keq = 0.25; Evalue = 39; % Next value here has concn units (M), then first order const % with s-1, then rest are dimensionless KMS = 0.1; % M units k2R = 1000; %s-1 SCRrat = 0.5; Erat = 0.5; SCR21rat = 0.02; SCSrat = 1; k24rat = 1; k2rat = 1; k24Srat = 1; k12Rrat = 2; k12Srat = 2; k34Rrat = 2; k34Srat = 2; if k2R < 0 % i.e. never, used to skip bits not wanted in current run dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, true); % For test end if k2R < 0 for S0 = [0.05 0.1 0.2 0.4 1] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end S0 = 0.2 end if k2R < 0 for Q0 = [0 0.25 0.5 1 2] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end Q0 = 0.5 end if k2R < 0 for Keq = [0.125 0.25 0.5 1 2 5 10 20 50 100] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end Keq = 0.25 %end %if k2R < 0 % never for Evalue = [9 19 39 100] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end Evalue = 39 %end %if k2R < 0 for KMS = [0.001 0.01 0.05 0.1 0.2 1 10] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end KMS = 0.1 end for k2R = [500 1000 2000] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end k2R = 1000 %end %if k2R < 0 for SCRrat = [0.25 0.5 1 2] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end SCRrat = 0.5 for Erat = [0.25 0.5 1 2] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end Erat = 0.5 for SCR21rat = [0.01 0.02 0.04 0.1 0.2 0.5 1 2] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end SCR21rat = 0.02 for SCSrat = [0.5 1 2] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end SCSrat = 1 for k24rat = [0.5 1 2] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end k24rat = 1 for k2rat = [0.5 1 2] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end k2rat = 1 for k24Srat = [0.5 1 2] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end k24Srat = 1 %end %if k2R < 0 for k12Rrat = [1 2 4] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end k12Rrat = 2 for k12Srat = [1 2 4] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end k12Srat = 2 for k34Rrat = [1 2 4] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end k34Rrat = 2 for k34Srat = [1 2 4] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end k34Srat = 2 %end if k2R<0 for Keq = [1 2 5 10 20] for S0 = [0.05 0.1 0.2 0.5 1 2] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end end Keq = 1 S0 = 0.2 %end for SCRrat = [0.1 0.2 0.5 1 2 5 10] %[0.1 0.2 0.5 1 2 5 10] for SCR21rat = [0.2 0.5 1 2 5 10] % 1 2 5 10] %[0.0005 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end end SCRrat = 0.5 SCR21rat = 1 for Erat = [0.1 0.2 0.5 1 2 5 10] %[0.1 0.2 0.5 1 2 5 10] for SCSrat = [0.1 0.2 0.5 1 2 5 10] % 1 2 5 10] %[0.0005 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5] dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false); end end Erat = 0.5 SCSrat = 1 end % Now the function, which doesn't actually return anything useful, just % does some jobs function dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, showplot) initrate = 0.005; %sets time scale, used to calculate EobySo % Now calculate rate constants, at first all on dimensioned basis k4R = k2R*k24rat; k2S = k2R/(k2rat*Evalue); k4S = k24Srat*k2S; denom1 = B0*KMS*SCR21rat*(k2rat + 1 + k2rat/k24Srat + 1/k24rat); SCRf = k2R*(B0*SCR21rat - Erat*KMS - KMS)/denom1; if SCRf < 0 error('SCRf negative') end SCRb = SCRf/SCRrat; SCSf = SCRf/Evalue; SCR2f = SCR21rat*SCRf; SCS2f = SCR2f/(Evalue*Erat); SCR2b = SCRf*SCR2f/(SCRb*Keq); SCSb = SCSf/SCSrat; SCS2b = SCRb*SCR2b*SCSf*SCS2f/(SCSb*SCRf*SCR2f); km1R = k12Rrat*k2R; k1R = SCRf*(1+k12Rrat); km2R = SCRb*(1+1/k12Rrat); km1S = k12Srat*k2S; k1S = SCSf*(1 + k12Srat); km2S = SCSb*(1 + 1/k12Srat); km3R = k34Rrat*k4R; k3R = SCR2f*(1 + k34Rrat); km4R = SCR2b*(1 + 1/k34Rrat); km3S = k34Srat*k4S; k3S = SCS2f*(1 + k34Srat); km4S = SCS2b*(1 + 1/k34Srat); EobySo = initrate*(1/(S0*SCRf)+1/k4R+1/k2R+1/(k4S*Evalue)+1/(k2S*Evalue)+ (1+Erat)/(B0*SCR21rat*SCRf)) if EobySo > 0.01 toobig = true; %error ('EobySo > 0.01') else toobig = false; end % Remove dimensions k1Rn = S0*k1R; k1Sn = S0*k1S; km2Rn = S0*km2R; km2Sn = S0*km2S; k3Rn = S0*k3R; k3Sn = S0*k3S; km4Rn = S0*km4R; km4Sn = S0*km4S; Sn0 = 1; % By definition of non-dimensional basis En0 = 1; % Same Bn0 = B0/S0; Qn0 = Q0/S0; % All other state variables start at 0 % Set up the values of state variables and simulation parameters % initial condition of state variables % En ESRn ESSn EsRn EsSn EPRn EPSn Sn Qn Bn PRn PSn Xinit = [En0 0 0 0 0 0 0 Sn0 Qn0 Bn0 0 0]'; % values of all rate constants used in integrating model k = [k1Rn k1Sn k2R k2S k3Rn k3Sn k4R k4S km1R km1S km2Rn km2Sn km3R km3S km4Rn km4Sn]'; % set the simulation time and sampling interval %tspan = 0:0.001:0.1; % start with transient %tspan = 0:0.4:20; % initial rate period tspan = 0:20:1000; % further progress %tspan = 0:20:3000; % further progress if ~toobig % otherwise skip this trial and wait for another % With normalised non-dimensional values for all variables and parameters, % probably similar AbsTol for all, so try keeping default 1e-6. AbsTol_value = 1e-6; options = odeset('AbsTol',AbsTol_value); % Run the ODE simulation %[t, x_s] = feval(@ode15s, @ODE_Prochiral, tspan, xs_init, options, k, EobySo); % This doesn't work any more when have extra parameter EobySo, but managed % to avoid feval call completely after viewing ode15s help %[t,y] = ode15s(@(t,y) odefcn(t,y,A,B), tspan, y0); Example from help [t, xout] = ode15s(@(t,x) ODE_ProchiralDiacid(t,x,k,EobySo), tspan, Xinit, options); % Help gives "[t,y] = ode15s(odefun,tspan,y0). Each row in the solution array %y corresponds to a value returned in column vector t." %disp(xout) %Yes, xout has a column for each state variable, and a row for each time, so x = transpose(xout); % x has one row for each state variable, one col for each time PRn = x(11,:); PSn = x(12,:); diff = PRn - PSn; conversion = PRn + PSn; ees = diff./conversion; else % is toobig showplot = false; t = '*'; conversion = 'EobySo too big'; ees = 'EobySo too big'; end % Prepare for output names = ["S0", "Q0", "B0", "Keq", "Eval", "KMS", "k2R", "SCRrat", "Erat", ... "SCR21rat", "SCSrat", "k24rat", "k2rat", "k24Srat", "k12Rrat", "k12Srat", ... "k34Rrat", "k34Srat"]; values = [S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ... SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat]; nonames = length(names); if length(values) ~= nonames error ('Different length names and values') end if showplot nolines = fix((nonames-1)/5) + 1; for lineno = 1:nolines firstctr = 5*(lineno-1) + 1; lastctr = 5*(lineno-1)+ 5; if lastctr > nonames lastctr = nonames; end lbl(lineno)=""; for ctr = firstctr:lastctr lbl(lineno) = strcat(lbl(lineno), names(ctr), "=", num2str(values(ctr)), " "); %thislbl = strcat(lbl(lineno), names(ctr), "=", num2str(values(ctr)), " ") end end figure (1) %plot(t,x(1,:),'b--',t,x(2,:),'r-',t,x(3,:),'g-.',t,x(4,:),'m-',t,x(5,:),'c-.','LineWidth',2.5); %plot(t,x(6,:),'b--',t,x(7,:),'r-',t,x(8,:),'g-.',t,x(9,:),'m-',t,x(10,:),'c-.','LineWidth',2.5); plot(t,conversion,'b--',t,ees,'r-',t,x(1,:),'g-.','LineWidth',2.5); hold on; %plot(t(b),a,'g.','MarkerSize',24); %text(t(b),a+0.05,max_text);xlabel('\fontsize{16}time/s'); % From help %dim = [.2 .5 .3 .3]; %annotation('textbox',dim,'String',str,'FitBoxToText','on'); % By trial and error, x and y run from origin in bottom left for ctr = 1:nolines dim = [0.15 0.67-0.07*ctr .3 .3]; annotation('textbox',dim,'String',lbl(ctr), 'FitBoxToText', 'on'); end ylabel('\fontsize{16}non-dimensional concentrations'); %title('\fontsize{14}time profile of state variables'); %legend('En','EDn','EtRn','EtSn','EQn'); %legend('Dn','Qn','Sn','PRn','PSn'); legend('Conv','ee_P','En'); end % Output_progress filename = 'C:\Users\ph\My Documents\A Edited\ProgressCurves.xlsx'; %rowno = int8(5) startcol = 5; firstdatacol = startcol + nonames + 2; rowno = FindBlankRow(filename, 5, firstdatacol); for ctr = 1:nonames writematrix(names(ctr), filename, 'Sheet', 1, 'Range', MakeCellRef(rowno, startcol+ctr-1)); writematrix(values(ctr), filename, 'Sheet', 1, 'Range', MakeCellRef(rowno+ 1, startcol+ctr-1)); end writematrix('time',filename, 'Sheet', 1, 'Range', MakeCellRef(rowno, firstdatacol)); writematrix('conv',filename, 'Sheet', 1, 'Range', MakeCellRef(rowno+1, firstdatacol)); writematrix('ee_P',filename, 'Sheet', 1, 'Range', MakeCellRef(rowno+2, firstdatacol)); %nextref = MakeCellRef(rowno,startcol+15); nextref = MakeCellRef(rowno,firstdatacol+1); writematrix(transpose(t),filename,'Sheet',1,'Range',nextref); nextref = MakeCellRef(rowno+1, firstdatacol+1); writematrix(conversion,filename,'Sheet',1,'Range',nextref); nextref = MakeCellRef(rowno+2, firstdatacol+1); writematrix(ees,filename,'Sheet',1,'Range',nextref); dummy = 1; %to prevent warning end function rowno = FindBlankRow(filename, startrow, colno) cellref = MakeCellRef(int8(startrow), int8(colno)); entries = readmatrix(filename, 'Sheet', 1, 'Range', strcat(cellref, ':', cellref)); norows = length(entries); rowno = startrow + norows; end function cellref = MakeCellRef(rowno, colno) firstletno = idivide(colno - 1, int8(26)); secondletno = mod(colno - 1, int8(26)) + 1; if firstletno > 0 cellref = char(firstletno + 64); else cellref = ''; end cellref = strcat(cellref,char(secondletno+64)); cellref = strcat(cellref,num2str(rowno)); end