Purpose: prepare ODEFUN for use with built-in Matlab solvers.
% Example: solve model using ode23
[ODEFUN,IC,elem_names]=dsDynasim2odefun(model);
[ODEFUN,IC,elem_names]=dsDynasim2odefun(dsPropagateParameters(dsPropagateFunctions(model)));
options=odeset('RelTol',1e-2,'AbsTol',1e-4,'InitialStep',.01);
[t,y]=ode23(ODEFUN,[0 100],IC,options);
figure; plot(t,y); legend(elem_names{:},'Location','EastOutside');
% Example: solve model manually using Euler method:
[ODEFUN,IC,elem_names]=dsDynasim2odefun(model);
dt=.01; t=0:dt:100;
y=zeros(length(t),length(IC));
y(1,:)=IC;
for i=2:length(t)
y(i,:)=y(i-1,:)+dt*ODEFUN(t,y(i-1,:)')';
end
figure; plot(t,y); legend(elem_names{:},'Location','EastOutside');
% without transposition:
dt=.01; t=0:dt:100;
y=zeros(length(IC),length(t));
y(:,1)=IC;
for i=2:length(t)
y(:,i)=y(:,i-1)+dt*ODEFUN(t,y(:,i-1));
end
figure; plot(t,y); legend(elem_names{:},'Location','EastOutside');
Note on implementation:
built-in solvers require ODEFUN to return state vector rows (i.e., cells
along rows), and they output state vectors with cells along columns
([y]=time x cells). In contrast, DynaSim ODEs/functions and output state
vectors have cells along columns. Therefore, for DynaSim models to be
compatible with built-in solvers, all state vectors must be transposed in
ODEFUN. This slows down simulation but cannot be avoided easily.
Author: Jason Sherfey, PhD <jssherfey@gmail.com>
Copyright (C) 2016 Jason Sherfey, Boston University, USA0001 function [ODEFUN,IC,elem_names] = dsDynasim2odefun(model, varargin) 0002 % Purpose: prepare ODEFUN for use with built-in Matlab solvers. 0003 % 0004 % % Example: solve model using ode23 0005 % [ODEFUN,IC,elem_names]=dsDynasim2odefun(model); 0006 % [ODEFUN,IC,elem_names]=dsDynasim2odefun(dsPropagateParameters(dsPropagateFunctions(model))); 0007 % options=odeset('RelTol',1e-2,'AbsTol',1e-4,'InitialStep',.01); 0008 % [t,y]=ode23(ODEFUN,[0 100],IC,options); 0009 % figure; plot(t,y); legend(elem_names{:},'Location','EastOutside'); 0010 % 0011 % % Example: solve model manually using Euler method: 0012 % [ODEFUN,IC,elem_names]=dsDynasim2odefun(model); 0013 % dt=.01; t=0:dt:100; 0014 % y=zeros(length(t),length(IC)); 0015 % y(1,:)=IC; 0016 % for i=2:length(t) 0017 % y(i,:)=y(i-1,:)+dt*ODEFUN(t,y(i-1,:)')'; 0018 % end 0019 % figure; plot(t,y); legend(elem_names{:},'Location','EastOutside'); 0020 % 0021 % % without transposition: 0022 % dt=.01; t=0:dt:100; 0023 % y=zeros(length(IC),length(t)); 0024 % y(:,1)=IC; 0025 % for i=2:length(t) 0026 % y(:,i)=y(:,i-1)+dt*ODEFUN(t,y(:,i-1)); 0027 % end 0028 % figure; plot(t,y); legend(elem_names{:},'Location','EastOutside'); 0029 % 0030 % Note on implementation: 0031 % built-in solvers require ODEFUN to return state vector rows (i.e., cells 0032 % along rows), and they output state vectors with cells along columns 0033 % ([y]=time x cells). In contrast, DynaSim ODEs/functions and output state 0034 % vectors have cells along columns. Therefore, for DynaSim models to be 0035 % compatible with built-in solvers, all state vectors must be transposed in 0036 % ODEFUN. This slows down simulation but cannot be avoided easily. 0037 % 0038 % Author: Jason Sherfey, PhD <jssherfey@gmail.com> 0039 % Copyright (C) 2016 Jason Sherfey, Boston University, USA 0040 0041 % Approach: 0042 % 1. evaluate params -> fixed_vars -> funcs 0043 % 2. evaluate ICs to get (# elems) per state var 0044 % 3. prepare state vector X 0045 % 4. replace state vars in ODEs by X 0046 % 5. combine X ODEs into ODEFUN 0047 0048 options=dsCheckOptions(varargin,{... 0049 'odefun_output','func_handle',{'func_handle','anonymous_func_string','func_body'},... 0050 'auto_gen_test_data_flag',0,{0,1},... 0051 },false); 0052 0053 if options.auto_gen_test_data_flag 0054 varargs = varargin; 0055 varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0; 0056 varargs(end+1:end+2) = {'unit_test_flag',1}; 0057 argin = [{model}, varargs]; % specific to this function 0058 end 0059 0060 % evaluate params -> fixed_vars -> funcs 0061 types={'parameters','fixed_variables','functions'}; 0062 for p=1:length(types) 0063 type=types{p}; 0064 if ~isempty(model.(type)) 0065 fields=fieldnames(model.(type)); 0066 for i=1:length(fields) 0067 val=model.(type).(fields{i}); 0068 if ~ischar(val) 0069 val=toString(val,'compact'); 0070 end 0071 % evaluate 0072 eval(sprintf('%s = %s;',fields{i},val)); 0073 % evalin('caller',sprintf('%s = %s;',fields{i},val)); 0074 % assignin('caller',fields{i},val); 0075 end 0076 end 0077 end 0078 0079 % evaluate ICs to get (# elems) per state var and set up generic state var X 0080 num_vars=length(model.state_variables); 0081 num_elems=zeros(1,num_vars); 0082 old_vars=model.state_variables; 0083 new_vars=cell(1,num_vars); 0084 new_inds=cell(1,num_vars); 0085 all_ICs=cell(1,num_vars); 0086 IC_names={}; 0087 state_var_index=0; 0088 0089 for i=1:num_vars 0090 var=model.state_variables{i}; 0091 0092 % check ICs for use of inital state_var value and put in proper starting value 0093 if regexp(model.ICs.(var), '_last') 0094 stateVars = regexp(model.ICs.(var), '([\w_]+)_last', 'tokens'); 0095 model.ICs.(var) = regexprep(model.ICs.(var), '_last', ''); 0096 0097 for iSVar = 1:length(stateVars) 0098 thisSvar = stateVars{iSVar}{1}; 0099 model.ICs.(var) = regexprep(model.ICs.(var), thisSvar, model.ICs.(thisSvar)); 0100 end 0101 end 0102 0103 % evaluate ICs to get (# elems) per state var 0104 ic=eval([model.ICs.(var) ';']); 0105 num_elems(i)=length(ic); 0106 0107 % set state var indices a variables for generic state vector X 0108 all_ICs{i}=ic; 0109 IC_names{i}=repmat({var},[1 num_elems(i)]); 0110 new_inds{i}=state_var_index+(1:length(ic)); 0111 new_vars{i}=sprintf('X(%g:%g)''',new_inds{i}(1),new_inds{i}(end)); 0112 state_var_index=state_var_index+length(ic); 0113 end 0114 0115 % prepare ODE system (comma-separated ODEs) 0116 ODEs=strtrim(struct2cell(model.ODEs)); 0117 idx=cellfun(@isempty,regexp(ODEs,';$')); % lines that need semicolons 0118 ODEs(idx)=cellfun(@(x)[x '&'],ODEs(idx),'uni',0); 0119 ODEs=[ODEs{:}]; % concatenate ODEs into a single string 0120 0121 % substitute in generic state vector X 0122 for i=1:num_vars 0123 ODEs=dsStrrep(ODEs,old_vars{i},new_vars{i}, '','', varargin{:}); 0124 end 0125 0126 % prepare outputs (function handle string, ICs, and element names for 0127 % mapping each X(i) to a particular state variable): 0128 elem_names=cat(2,IC_names{:}); 0129 0130 switch options.odefun_output 0131 case 'func_handle' 0132 ODEs=strrep(ODEs,'&',','); % replace semicolons by commas 0133 ODEs(end) = []; %remove trailing comma 0134 ODEFUN = eval(['@(t,X) [' ODEs ']'';']); 0135 case 'anonymous_func_string' 0136 ODEs=strrep(ODEs,'&',','); % replace semicolons by commas 0137 ODEs(end) = []; %remove trailing comma 0138 ODEFUN = ['@(t,X) [' ODEs ']'';']; 0139 case 'func_body' 0140 ODEs=strrep(ODEs,'&',',...\n'); % replace semicolons by commas with newline 0141 ODEFUN = ODEs; 0142 end 0143 0144 IC=cat(2,all_ICs{:})'; 0145 0146 %% auto_gen_test_data_flag argout 0147 if options.auto_gen_test_data_flag 0148 argout = {ODEFUN, IC, elem_names}; % specific to this function 0149 0150 dsUnitSaveAutoGenTestData(argin, argout); 0151 end