Generated Code
The following is matlab code generated by the CellML API from this CellML file. (Back to language selection)
The raw code is available.
function [VOI, STATES, ALGEBRAIC, CONSTANTS] = mainFunction() % This is the "main function". In Matlab, things work best if you rename this function to match the filename. [VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel(); end function [algebraicVariableCount] = getAlgebraicVariableCount() % Used later when setting a global variable with the number of algebraic variables. % Note: This is not the "main method". algebraicVariableCount =15; end % There are a total of 3 entries in each of the rate and state variable arrays. % There are a total of 21 entries in the constant variable array. % function [VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel() % Create ALGEBRAIC of correct size global algebraicVariableCount; algebraicVariableCount = getAlgebraicVariableCount(); % Initialise constants and state variables [INIT_STATES, CONSTANTS] = initConsts; % Set timespan to solve over tspan = [0, 10]; % Set numerical accuracy options for ODE solver options = odeset('RelTol', 1e-06, 'AbsTol', 1e-06, 'MaxStep', 1); % Solve model with ODE solver [VOI, STATES] = ode15s(@(VOI, STATES)computeRates(VOI, STATES, CONSTANTS), tspan, INIT_STATES, options); % Compute algebraic variables [RATES, ALGEBRAIC] = computeRates(VOI, STATES, CONSTANTS); ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, VOI); % Plot state variables against variable of integration [LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends(); figure(); plot(VOI, STATES); xlabel(LEGEND_VOI); l = legend(LEGEND_STATES); set(l,'Interpreter','none'); end function [LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends() LEGEND_STATES = ''; LEGEND_ALGEBRAIC = ''; LEGEND_VOI = ''; LEGEND_CONSTANTS = ''; LEGEND_VOI = strpad('t in component environment (second)'); LEGEND_ALGEBRAIC(:,8) = strpad('Pi in component TempCDa (UnitP)'); LEGEND_ALGEBRAIC(:,12) = strpad('Pi in component TempR (UnitP)'); LEGEND_ALGEBRAIC(:,11) = strpad('Qo in component TempRC (UnitQ)'); LEGEND_ALGEBRAIC(:,13) = strpad('Qo in component TempCDv (UnitQ)'); LEGEND_ALGEBRAIC(:,4) = strpad('Pi in component TempCDv (UnitP)'); LEGEND_ALGEBRAIC(:,10) = strpad('Qo in component TempCDa (UnitQ)'); LEGEND_CONSTANTS(:,1) = strpad('CVao in component ParaLeftHeart (UnitCV)'); LEGEND_ALGEBRAIC(:,3) = strpad('E in component EVentricle (UnitE)'); LEGEND_STATES(:,1) = strpad('V in component TempCDv (UnitV)'); LEGEND_CONSTANTS(:,2) = strpad('PlvIni in component ParaLeftHeart (UnitP)'); LEGEND_CONSTANTS(:,3) = strpad('VlvIni in component ParaLeftHeart (UnitV)'); LEGEND_ALGEBRAIC(:,14) = strpad('Tao in component TempCDv (dimensionless)'); LEGEND_CONSTANTS(:,4) = strpad('Vlv0 in component ParaLeftHeart (UnitV)'); LEGEND_CONSTANTS(:,5) = strpad('CVmi in component ParaLeftHeart (UnitCV)'); LEGEND_ALGEBRAIC(:,7) = strpad('E in component EAtrium (UnitE)'); LEGEND_STATES(:,2) = strpad('V in component TempCDa (UnitV)'); LEGEND_CONSTANTS(:,6) = strpad('PlaIni in component ParaLeftHeart (UnitP)'); LEGEND_CONSTANTS(:,7) = strpad('VlaIni in component ParaLeftHeart (UnitV)'); LEGEND_ALGEBRAIC(:,9) = strpad('Tao in component TempCDa (dimensionless)'); LEGEND_CONSTANTS(:,8) = strpad('Vla0 in component ParaLeftHeart (UnitV)'); LEGEND_CONSTANTS(:,9) = strpad('ElvMax in component ParaLeftHeart (UnitE)'); LEGEND_CONSTANTS(:,10) = strpad('ElvMin in component ParaLeftHeart (UnitE)'); LEGEND_CONSTANTS(:,11) = strpad('T in component ParaLeftHeart (second)'); LEGEND_CONSTANTS(:,12) = strpad('Ts1 in component ParaLeftHeart (dimensionless)'); LEGEND_CONSTANTS(:,13) = strpad('Ts2 in component ParaLeftHeart (dimensionless)'); LEGEND_ALGEBRAIC(:,1) = strpad('mt in component EVentricle (second)'); LEGEND_ALGEBRAIC(:,2) = strpad('et in component EVentricle (dimensionless)'); LEGEND_CONSTANTS(:,14) = strpad('ElaMax in component ParaLeftHeart (UnitE)'); LEGEND_CONSTANTS(:,15) = strpad('ElaMin in component ParaLeftHeart (UnitE)'); LEGEND_CONSTANTS(:,16) = strpad('Tpwb in component ParaLeftHeart (dimensionless)'); LEGEND_CONSTANTS(:,17) = strpad('Tpww in component ParaLeftHeart (dimensionless)'); LEGEND_ALGEBRAIC(:,5) = strpad('mt in component EAtrium (second)'); LEGEND_ALGEBRAIC(:,6) = strpad('et in component EAtrium (dimensionless)'); LEGEND_STATES(:,3) = strpad('Pi in component TempRC (UnitP)'); LEGEND_ALGEBRAIC(:,15) = strpad('Qo in component TempR (UnitQ)'); LEGEND_CONSTANTS(:,18) = strpad('Rc in component ParaWestkessel (UnitR)'); LEGEND_CONSTANTS(:,19) = strpad('Rv in component ParaWestkessel (UnitR)'); LEGEND_CONSTANTS(:,20) = strpad('Cv in component ParaWestkessel (UnitC)'); LEGEND_CONSTANTS(:,21) = strpad('P0v in component ParaWestkessel (UnitP)'); LEGEND_RATES(:,1) = strpad('d/dt V in component TempCDv (UnitV)'); LEGEND_RATES(:,2) = strpad('d/dt V in component TempCDa (UnitV)'); LEGEND_RATES(:,3) = strpad('d/dt Pi in component TempRC (UnitP)'); LEGEND_STATES = LEGEND_STATES'; LEGEND_ALGEBRAIC = LEGEND_ALGEBRAIC'; LEGEND_RATES = LEGEND_RATES'; LEGEND_CONSTANTS = LEGEND_CONSTANTS'; end function [STATES, CONSTANTS] = initConsts() VOI = 0; CONSTANTS = []; STATES = []; ALGEBRAIC = []; CONSTANTS(:,1) = 350.; CONSTANTS(:,2) = 1.0; CONSTANTS(:,3) = 5.0; CONSTANTS(:,4) = 500; CONSTANTS(:,5) = 400.; CONSTANTS(:,6) = 1.0; CONSTANTS(:,7) = 4.0; CONSTANTS(:,8) = 20; CONSTANTS(:,9) = 2.5; CONSTANTS(:,10) = 0.1; CONSTANTS(:,11) = 1.0; CONSTANTS(:,12) = 0.3; CONSTANTS(:,13) = 0.45; CONSTANTS(:,14) = 0.25; CONSTANTS(:,15) = 0.15; CONSTANTS(:,16) = 0.92; CONSTANTS(:,17) = 0.09; CONSTANTS(:,18) = 0.005; CONSTANTS(:,19) = 0.65; CONSTANTS(:,20) = 2.6; CONSTANTS(:,21) = 0.; STATES(:,1) = CONSTANTS(:,4); STATES(:,2) = CONSTANTS(:,8); STATES(:,3) = CONSTANTS(:,21); if (isempty(STATES)), warning('Initial values for states not set');, end end function [RATES, ALGEBRAIC] = computeRates(VOI, STATES, CONSTANTS) global algebraicVariableCount; statesSize = size(STATES); statesColumnCount = statesSize(2); if ( statesColumnCount == 1) STATES = STATES'; ALGEBRAIC = zeros(1, algebraicVariableCount); utilOnes = 1; else statesRowCount = statesSize(1); ALGEBRAIC = zeros(statesRowCount, algebraicVariableCount); RATES = zeros(statesRowCount, statesColumnCount); utilOnes = ones(statesRowCount, 1); end ALGEBRAIC(:,5) = VOI - CONSTANTS(:,11).*floor(VOI./CONSTANTS(:,11)); ALGEBRAIC(:,6) = piecewise({ALGEBRAIC(:,5)>=0.00000&ALGEBRAIC(:,5)<= ((CONSTANTS(:,16)+CONSTANTS(:,17)) - 1.00000).*CONSTANTS(:,11), 1.00000 - cos(( 2.00000.*3.14159.*(ALGEBRAIC(:,5) - (CONSTANTS(:,16) - 1.00000).*CONSTANTS(:,11)))./( CONSTANTS(:,17).*CONSTANTS(:,11))) , ALGEBRAIC(:,5)> ((CONSTANTS(:,16)+CONSTANTS(:,17)) - 1.00000).*CONSTANTS(:,11)&ALGEBRAIC(:,5)<= CONSTANTS(:,16).*CONSTANTS(:,11), 0.00000 , ALGEBRAIC(:,5)> CONSTANTS(:,16).*CONSTANTS(:,11)&ALGEBRAIC(:,5)<=CONSTANTS(:,11), 1.00000 - cos(( 2.00000.*3.14159.*(ALGEBRAIC(:,5) - CONSTANTS(:,16).*CONSTANTS(:,11)))./( CONSTANTS(:,17).*CONSTANTS(:,11))) }, NaN); ALGEBRAIC(:,7) = CONSTANTS(:,15)+( ALGEBRAIC(:,6).*(CONSTANTS(:,14) - CONSTANTS(:,15)))./2.00000; ALGEBRAIC(:,8) = CONSTANTS(:,6)+ ALGEBRAIC(:,7).*(STATES(:,2) - CONSTANTS(:,7)); ALGEBRAIC(:,11) = (STATES(:,3) - ALGEBRAIC(:,8))./CONSTANTS(:,19); ALGEBRAIC(:,1) = VOI - CONSTANTS(:,11).*floor(VOI./CONSTANTS(:,11)); ALGEBRAIC(:,2) = piecewise({ALGEBRAIC(:,1)>=0.00000&ALGEBRAIC(:,1)<= CONSTANTS(:,12).*CONSTANTS(:,11), 1.00000 - cos(( 3.14159.*ALGEBRAIC(:,1))./( CONSTANTS(:,12).*CONSTANTS(:,11))) , ALGEBRAIC(:,1)> CONSTANTS(:,12).*CONSTANTS(:,11)&ALGEBRAIC(:,1)<= CONSTANTS(:,13).*CONSTANTS(:,11), 1.00000+cos(( 3.14159.*(ALGEBRAIC(:,1) - CONSTANTS(:,12).*CONSTANTS(:,11)))./( (CONSTANTS(:,13) - CONSTANTS(:,12)).*CONSTANTS(:,11))) , ALGEBRAIC(:,1)> CONSTANTS(:,13).*CONSTANTS(:,11)&ALGEBRAIC(:,1)<CONSTANTS(:,11), 0.00000 }, NaN); ALGEBRAIC(:,3) = CONSTANTS(:,10)+( ALGEBRAIC(:,2).*(CONSTANTS(:,9) - CONSTANTS(:,10)))./2.00000; ALGEBRAIC(:,4) = CONSTANTS(:,2)+ ALGEBRAIC(:,3).*(STATES(:,1) - CONSTANTS(:,3)); ALGEBRAIC(:,9) = piecewise({ALGEBRAIC(:,8)>=ALGEBRAIC(:,4), 1.00000 , ALGEBRAIC(:,8)<ALGEBRAIC(:,4), 0.00000 }, NaN); ALGEBRAIC(:,10) = piecewise({ALGEBRAIC(:,8)>=ALGEBRAIC(:,4), CONSTANTS(:,5).*ALGEBRAIC(:,9).*power(abs(ALGEBRAIC(:,8) - ALGEBRAIC(:,4)), 0.500000) , ALGEBRAIC(:,8)<ALGEBRAIC(:,4), -1.00000.*CONSTANTS(:,5).*ALGEBRAIC(:,9).*power(abs(ALGEBRAIC(:,4) - ALGEBRAIC(:,8)), 0.500000) }, NaN); RATES(:,2) = ALGEBRAIC(:,11) - ALGEBRAIC(:,10); [CONSTANTS, STATES, ALGEBRAIC] = rootfind_0(VOI, CONSTANTS, STATES, ALGEBRAIC); RATES(:,1) = ALGEBRAIC(:,10) - ALGEBRAIC(:,13); ALGEBRAIC(:,15) = ALGEBRAIC(:,13); RATES(:,3) = (ALGEBRAIC(:,15) - ALGEBRAIC(:,11))./CONSTANTS(:,20); RATES = RATES'; end % Calculate algebraic variables function ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, VOI) statesSize = size(STATES); statesColumnCount = statesSize(2); if ( statesColumnCount == 1) STATES = STATES'; utilOnes = 1; else statesRowCount = statesSize(1); utilOnes = ones(statesRowCount, 1); end ALGEBRAIC(:,5) = VOI - CONSTANTS(:,11).*floor(VOI./CONSTANTS(:,11)); ALGEBRAIC(:,6) = piecewise({ALGEBRAIC(:,5)>=0.00000&ALGEBRAIC(:,5)<= ((CONSTANTS(:,16)+CONSTANTS(:,17)) - 1.00000).*CONSTANTS(:,11), 1.00000 - cos(( 2.00000.*3.14159.*(ALGEBRAIC(:,5) - (CONSTANTS(:,16) - 1.00000).*CONSTANTS(:,11)))./( CONSTANTS(:,17).*CONSTANTS(:,11))) , ALGEBRAIC(:,5)> ((CONSTANTS(:,16)+CONSTANTS(:,17)) - 1.00000).*CONSTANTS(:,11)&ALGEBRAIC(:,5)<= CONSTANTS(:,16).*CONSTANTS(:,11), 0.00000 , ALGEBRAIC(:,5)> CONSTANTS(:,16).*CONSTANTS(:,11)&ALGEBRAIC(:,5)<=CONSTANTS(:,11), 1.00000 - cos(( 2.00000.*3.14159.*(ALGEBRAIC(:,5) - CONSTANTS(:,16).*CONSTANTS(:,11)))./( CONSTANTS(:,17).*CONSTANTS(:,11))) }, NaN); ALGEBRAIC(:,7) = CONSTANTS(:,15)+( ALGEBRAIC(:,6).*(CONSTANTS(:,14) - CONSTANTS(:,15)))./2.00000; ALGEBRAIC(:,8) = CONSTANTS(:,6)+ ALGEBRAIC(:,7).*(STATES(:,2) - CONSTANTS(:,7)); ALGEBRAIC(:,11) = (STATES(:,3) - ALGEBRAIC(:,8))./CONSTANTS(:,19); ALGEBRAIC(:,1) = VOI - CONSTANTS(:,11).*floor(VOI./CONSTANTS(:,11)); ALGEBRAIC(:,2) = piecewise({ALGEBRAIC(:,1)>=0.00000&ALGEBRAIC(:,1)<= CONSTANTS(:,12).*CONSTANTS(:,11), 1.00000 - cos(( 3.14159.*ALGEBRAIC(:,1))./( CONSTANTS(:,12).*CONSTANTS(:,11))) , ALGEBRAIC(:,1)> CONSTANTS(:,12).*CONSTANTS(:,11)&ALGEBRAIC(:,1)<= CONSTANTS(:,13).*CONSTANTS(:,11), 1.00000+cos(( 3.14159.*(ALGEBRAIC(:,1) - CONSTANTS(:,12).*CONSTANTS(:,11)))./( (CONSTANTS(:,13) - CONSTANTS(:,12)).*CONSTANTS(:,11))) , ALGEBRAIC(:,1)> CONSTANTS(:,13).*CONSTANTS(:,11)&ALGEBRAIC(:,1)<CONSTANTS(:,11), 0.00000 }, NaN); ALGEBRAIC(:,3) = CONSTANTS(:,10)+( ALGEBRAIC(:,2).*(CONSTANTS(:,9) - CONSTANTS(:,10)))./2.00000; ALGEBRAIC(:,4) = CONSTANTS(:,2)+ ALGEBRAIC(:,3).*(STATES(:,1) - CONSTANTS(:,3)); ALGEBRAIC(:,9) = piecewise({ALGEBRAIC(:,8)>=ALGEBRAIC(:,4), 1.00000 , ALGEBRAIC(:,8)<ALGEBRAIC(:,4), 0.00000 }, NaN); ALGEBRAIC(:,10) = piecewise({ALGEBRAIC(:,8)>=ALGEBRAIC(:,4), CONSTANTS(:,5).*ALGEBRAIC(:,9).*power(abs(ALGEBRAIC(:,8) - ALGEBRAIC(:,4)), 0.500000) , ALGEBRAIC(:,8)<ALGEBRAIC(:,4), -1.00000.*CONSTANTS(:,5).*ALGEBRAIC(:,9).*power(abs(ALGEBRAIC(:,4) - ALGEBRAIC(:,8)), 0.500000) }, NaN); ALGEBRAIC(:,15) = ALGEBRAIC(:,13); end % Functions required for solving differential algebraic equation function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_0(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN) ALGEBRAIC = ALGEBRAIC_IN; CONSTANTS = CONSTANTS_IN; STATES = STATES_IN; global initialGuess_0; if (length(initialGuess_0) ~= 3), initialGuess_0 = [0.1,0.1,0.1];, end options = optimset('Display', 'off', 'TolX', 1E-6); if length(VOI) == 1 residualfn = @(algebraicCandidate)residualSN_0(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES); soln = fsolve(residualfn, initialGuess_0, options); initialGuess_0 = soln; ALGEBRAIC(:,12) = soln(1); ALGEBRAIC(:,13) = soln(2); ALGEBRAIC(:,14) = soln(3); else SET_ALGEBRAIC(:,12) = logical(1); SET_ALGEBRAIC(:,13) = logical(1); SET_ALGEBRAIC(:,14) = logical(1); for i=1:length(VOI) residualfn = @(algebraicCandidate)residualSN_0(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:)); soln = fsolve(residualfn, initialGuess_0, options); initialGuess_0 = soln; TEMP_ALGEBRAIC(:,12) = soln(1); TEMP_ALGEBRAIC(:,13) = soln(2); TEMP_ALGEBRAIC(:,14) = soln(3); ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC); end end end function resid = residualSN_0(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES) ALGEBRAIC(:,12) = algebraicCandidate(1); ALGEBRAIC(:,13) = algebraicCandidate(2); ALGEBRAIC(:,14) = algebraicCandidate(3); resid(1) = ALGEBRAIC(:,14) - piecewise({ALGEBRAIC(:,4)>=ALGEBRAIC(:,12), 1.00000 , ALGEBRAIC(:,4)<ALGEBRAIC(:,12), 0.00000 }, NaN); resid(2) = ALGEBRAIC(:,13) - piecewise({ALGEBRAIC(:,4)>=ALGEBRAIC(:,12), CONSTANTS(:,1).*ALGEBRAIC(:,14).*power(abs(ALGEBRAIC(:,4) - ALGEBRAIC(:,12)), 0.500000) , ALGEBRAIC(:,4)<ALGEBRAIC(:,12), -1.00000.*CONSTANTS(:,1).*ALGEBRAIC(:,14).*power(abs(ALGEBRAIC(:,12) - ALGEBRAIC(:,4)), 0.500000) }, NaN); resid(3) = ALGEBRAIC(:,12) - (STATES(:,3)+ CONSTANTS(:,18).*ALGEBRAIC(:,13)); end % Compute result of a piecewise function function x = piecewise(cases, default) set = [0]; for i = 1:2:length(cases) if (length(cases{i+1}) == 1) x(cases{i} & ~set,:) = cases{i+1}; else x(cases{i} & ~set,:) = cases{i+1}(cases{i} & ~set); end set = set | cases{i}; if(set), break, end end if (length(default) == 1) x(~set,:) = default; else x(~set,:) = default(~set); end end % Pad out or shorten strings to a set length function strout = strpad(strin) req_length = 160; insize = size(strin,2); if insize > req_length strout = strin(1:req_length); else strout = [strin, blanks(req_length - insize)]; end end