Can I use numeric::odesolve as a replacement to ode45?

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.

Answer 1

[Walter Roberson]

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.

[srikumar C Gopalakrishnan]

I want to implement function in ODESolve with values specified as input vector and initial conditions. Is there a way to implement the input in this manner ?. Any input on this is highly appreciated !!!

[Walter Roberson]

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

Answer 2

[srikumar C Gopalakrishnan]

function ActualImplementODEFcn
A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode45(@(t,y) odefcn(t,y,A,B), tspan, y0);
save tvalue t
save yvalue y
end

function dydt = odefcn(t,y,A,B)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = (A/B)*t.*y(1);
end

Here is a function which is similar to the one I am implementing. This code fragment is from MATLAB help to implement ode45 and odefcn calculates the derivatives which has to be used as input to calculate y. I am solving 2nd order differential equation, so I need 2 derivatives 1st derivative and 2nd derivative.

Answer 3

[srikumar C Gopalakrishnan]

function FunctionHandlesTaking
clear all
clc
syms f(tV)
syms g(tV)
syms dydtfirst
A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode45(@(t,y) odefcn(t,y,A,B), tspan, y0);
y=y0;
tempt=t;
yvalues=zeros(size(t,1),2);
yvalues(1,1)=0;
yvalues(1,2)=0.01;
for ii=1:1:length(tempt)
t=tempt(ii);
dydt=odefcn(t,y,A,B);
%eqns={diff(dydt,t)==dydt(2),dydt == dydt(1),y(1) == 0.01};
if (ii<length(tempt))
tiplus1 = tempt(ii+1);
else
tiplus1=max(tempt)+(max(tempt)/length(tempt));    
end
eqns={diff(f(tV),tV) == dydt(1),f(0)==y(1)};
limitvalue={t,tiplus1};
numdig = 16; 
digits(numdig);
symLIST = @(varargin)feval(symengine,'DOM_LIST',varargin{:});
symRANGE = @(a,b) feval(symengine, '_range', a, b);
IVP = symLIST(eqns{:});
fields = symLIST(f(tV));
IVP = feval(symengine,'numeric::ode2vectorfield', IVP, fields);
syms CK45
method = CK45;
y(1)=feval(symengine, 'numeric::odesolve', IVP(1), symRANGE(limitvalue{1,1}, limitvalue{1,2}), IVP(3), method);
yvalues(ii+1,1)=y(1);
eqns1={diff(g(tV),tV)== dydt(2), g(0)==y(2)};
numdig = 16;
digits(numdig);
symLIST = @(varargin)feval(symengine,'DOM_LIST',varargin{:});
symRANGE = @(c,d) feval(symengine, '_range', c, d);
IVP1 = symLIST(eqns1{:});
fields1 = symLIST(g(tV));
IVP1 = feval(symengine,'numeric::ode2vectorfield', IVP1, fields1);
syms CK45
method = CK45;
y(2)=feval(symengine, 'numeric::odesolve', IVP1(1), symRANGE(limitvalue{1,1}, limitvalue{1,2}), IVP1(3), method);
yvalues(ii+1,2)=y(2);
end
save yvalues yvalues
save tempt tempt
end

function dydt = odefcn(t,y,A,B)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = (A/B)*t.*y(1);
end

Above is my implementation of calling Cash Karp method from matlab. The code works fine please let me know if there are any issues in the manner I have implemented this CASH Karp calling from MATLAB. I compared the results generated by ODE45 Dormand Prince method and Cash Karp method, I am getting difference in the integrated final values, I am getting a maximum variation is 32%. I know the weightage used in Dormand Prince and Cash Karp methods are different, but I am not sure whether that will create so much of variation. I am not sure if I made any errors in the way we call Cash Karp Method from Symbolic Tool Box. Any comments and help on this implementation will be highly appreciated !.

[srikumar C Gopalakrishnan]

I have attached M-File which I am using with this answers. I am calculating both derivative and the variable using Cash-Karp Method. I am solving 2nd order differential equation using ODE Cash Karp method. How to implement variable time step in Cash-Karp method?

Answer 4

[srikumar C Gopalakrishnan]

[srikumar C Gopalakrishnan]

If we observer carefully the Cash Karp Y Prime Value matches with the ODE45 Y Value. I dont know why?

[Karan Gill]

Why do you need to use numeric::odesolve? Why isn't ode45 good enough?

[Walter Roberson]

User needs to specifically compare with Cash Karp which the numeric ode routines do not support.

[srikumar C Gopalakrishnan]

Cash-Karp method is said to handle a differential equation with rapidly varying right hand side. The Dormand Prince method (ODE45) which Iam using runs into stability problems. so I decided to compare the Dormand Prince method with Cash Karp method. This Cash Karp solver is not available in the numeric ode routines, so I have to use Symbolic tool box. I have compared Dormand Prince Solution with Cash Karp method and analytical method together.