https://github.com/Data2Dynamics/d2d
Tip revision: 83d6eedd5f6c33ae94dbce6c4e620d2f1bb57b1f authored by Joep Vanlier on 09 April 2015, 08:41:55 UTC
Forgot to update C version number
Forgot to update C version number
Tip revision: 83d6eed
arSolveSS.m
% solve steady state condition
%
% arSolveSS(m, fu, fextra)
%
% m: model index
% fu: custom initial input
% fextra: additional equations
function arSolveSS(m, fu, fextra)
global ar
fprintf('solving steady state condition...\n');
rhs = sym(ar.model(m).fx);
if(exist('fu','var'))
rhs = mysubs(rhs, ar.model(m).sym.u, sym(fu));
elseif(~isempty(ar.model(m).fu))
rhs = mysubs(rhs, ar.model(m).sym.u, sym(ar.model(m).fu));
end
rhs = mysubs(rhs, ar.model(m).sym.x, ar.model(m).sym.px0);
qnonzero = rhs~=0;
fstrings = cell(1, sum(logical(qnonzero)));
dxdt = cell(1, sum(logical(qnonzero)));
fcount = 1;
for j=find(qnonzero)'
fstrings{fcount} = sprintf('%s = 0', char(rhs(j)));
dxdt{fcount} = ar.model(m).x{j};
fcount = fcount + 1;
end
fp = cell(1, sum(logical(qnonzero)));
fcount = 1;
for j=find(qnonzero)'
fp{fcount} = ar.model(m).px0{j};
fcount = fcount + 1;
end
fprintf('\nsteady state equations (%i):\n', length(fstrings));
for j=1:length(fstrings)
disp(['d' dxdt{j} '/dt = ' fstrings{j}])
end
fprintf('\nsolving for (%i):\n', length(fp));
for j=1:length(fp)
disp(fp{j})
end
fprintf('\n---\n');
if(exist('fextra','var'))
ar.model(m).ss_solution = solve(fstrings{:}, fextra{:}, fp{:});
else
ar.model(m).ss_solution = solve(fstrings{:}, fp{:});
end
if(~isempty(ar.model(m).ss_solution))
fprintf('---\n');
tmpstr = getfield(ar.model(m).ss_solution, fp{1}); %#ok<GFLD>
for jj = 1:length(tmpstr)
fprintf('\n%s Solution %i:\n', ar.config.comment_string, jj);
for j = 1:length(fp)
tmpstr = getfield(ar.model(m).ss_solution, fp{j}); %#ok<GFLD>
%fprintf('%s\t "%s"\n', fp{j}, char(simple(tmpstr(jj))));
fprintf('%s\t "%s"\n', fp{j}, char(tmpstr(jj)));
end
end
end
% better subs
function out = mysubs(in, old, new)
if(~isnumeric(in) && ~isempty(old) && ~isempty(findsym(in)))
matVer = ver('MATLAB');
if(str2double(matVer.Version)>=8.1)
out = subs(in, old(:), new(:));
else
out = subs(in, old(:), new(:), 0);
end
else
out = in;
end