Can I use numeric::odesolve as a replacement to ode45?
7 views (last 30 days)
Show older comments
srikumar C Gopalakrishnan
on 26 Oct 2016
Commented: srikumar C Gopalakrishnan
on 23 Nov 2016
I am using ode45 to solve nonlinear differential equations. Can I use numeric::odesolve with CK45 as a substitute to ode45?. I did not find any material on using odesolve with Cash Karp as a replacement to ode45. Is there any differences in the inputs that has to be given to ode45 and CK45 in numeric odesolve?. Is there any thing which has to be kept in mind while converting the ode45 statement to CK45 numeric odesolve? Any inputs on this will be highly appreciated.
0 Comments
Accepted Answer
Walter Roberson
on 26 Oct 2016
It is possible, by using some magic. See for example,
syms y(t)
%Boundary conditions are required
eqns = {diff(y(t),t) == t*sin(y(t)), y(0) == 2}; %use a cell array!
numdig = 16; %number of decimal places to use
digits(numdig);
symLIST = @(varargin)feval(symengine,'DOM_LIST',varargin{:});
symRANGE = @(a,b) feval(symengine, '_range', a, b);
IVP = symLIST(eqns{:});
fields = symLIST(y(t));
IVP = feval(symengine,'numeric::ode2vectorfield', IVP, fields);
syms CK45
method = CK45;
feval(symengine, 'numeric::odesolve', IVP(1), symRANGE(IVP(2),10), IVP(3), method)
You can omit the method parameter from the call.
2 Comments
Walter Roberson
on 1 Nov 2016
It is important that you use the symLIST and symRANGE I show, or their equivalents. numeric::odesolve takes its parameters in MuPAD structures that you cannot normally create in MATLAB, and will not work if you pass in variables created in MATLAB. The printable representation might look exactly the same in some cases, but the structure is important.
For example in places you cannot just pass in 2 and you cannot just pass in sym(2), you need MuPAD's list form, which inside the MuPAD engine is designated as [2] but that is not the same as MATLAB's [2] or [sym(2)] even though MATLAB will display sym(2) and evalin(symengine,'[2]') the same way.
I found it easier not to call numeric::odesolve as the first step: constructing the function in the proper way was getting too hard unless I resorted to using evalin(symengine, sprintf(.....)) to build the appropriate string, in a way that was difficult to generalize to all functions. Calling numeric::ode2vectorfield to prepare the function was much easier when I found it.
You need to create a series of equations, probably including at least diff(), for numeric::odesolve to be used. I do not find any equations in what you posted here: I see an expression and I see some numeric values, but no equations. You can create a cell array of equations from symbols and values by arrangements such as... ummm...
bv = 1:8;
N = length(bv);
syms t
funs_cell = arrayfun( @(N) evalin(symengine, sprintf('f%d(t)', N)), 1:N, 'Uniform', 0);
funs_prime_cell = cellfun(@(F) diff(F,t), funs_cell, 'Uniform', 0);
rhs_cell = cellfun(@(F) F * sin(F), funs_cell, 'Uniform', 0);
eqns_cell = cellfun(@eq, funs_prime_cell, rhs_cell, 'Uniform', 0);
bc_cell = cellfun(@(F, BV) subs(F, t, 0) == BV, funs_cell, num2cell(bv), 'Uniform', 0);
eqns = [eqns_cell, bc_cell];
... and continue from the numdig point
These lines are constructing f1(t), f2(t) ..., and then diff(f1(t),t), diff(f2(t),t) ..., and then f1(t)*sin(f1(t)), f2(t)*sin(f2(t)) ..., then putting those together to get diff(f1(t),t)=f1(t)*sin(f1(t)), diff(f2(t),t)=f2(t)*sin(f2(t)) ... and then constructing f1(0) = bv(1), f2(0) = bv(2), ... and then putting all the diff() and f*(0) together .
You might well be tempted to just sprintf() the parts together, something like
for K = 1 : N
eqns_cell{K} = evalin(symengine, sprintf('diff(f%d(t),t)=f%d(t)*sin(f%d(t))', K, K, K) );
bc_cell{K} = evalin(symengine, sprintf('f%d(0)=%g', K, bv(K)));
end
eqns = [eqns_cell, bc_cell];
... and continue from the numdig point
More Answers (3)
srikumar C Gopalakrishnan
on 13 Nov 2016
4 Comments
Walter Roberson
on 21 Nov 2016
User needs to specifically compare with Cash Karp which the numeric ode routines do not support.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!