Variable precision arithmetic in ode45 and others?

Is it possible to use vpa in ordinary differential equation solvers like ode45?

Answer 1

[Walter Roberson]

Sort of. The internal symbolic engine has two numeric ode routines. At the moment I do not know if there is a matlab level interface to those routines.

You might need to use feval(symengine) to access the routines.

http://www.cfm.brown.edu/people/dobrush/am33/MuPad/MuPad22.html

[Walter Roberson]

I confirm that the functions continue to work in R2020a. For example,

S2 = "IVP2 := {y''(t) = t*y'(t), y(0) = 0, y'(0) = 12}; fields2 := [y(t), y'(t)];ODE2 := numeric::ode2vectorfield(IVP2, fields2); numApprox := numeric::odesolve2(ODE2); numApprox(4)";
evalin(symengine, S2)

ans =
[ 9672.4799720575046729200296634973, 35771.495844500739296923105193435]

I will poke around a bit and see what kind of interface I can make at the MATLAB level.

[Walter Roberson]

The below works so-far, so I am correctly constructing the IVP and fields. However, when I break S3 up into two separate commands, I get an error message after the first one... Ah, it appears to have to do with the MATLAB level trying to parse the return values, so I might be able to work with that.

Symlist = @(LIST) feval(symengine, 'DOM_LIST@op', LIST);
Symassign = @(VAR, VALUE) feval(symengine, '_assign', VAR, VALUE);

syms y(t)
Dy = diff(y);
D2y = diff(Dy);
IVP3 = Symassign('IVP3', Symlist([D2y == t*Dy, y(0)==0, Dy(0)==2]));
fields3 = Symassign('fields3', Symlist([y, Dy]));

S3 = "[ODE3, t0, Y0] := [numeric::ode2vectorfield(IVP3, fields3)]; numeric::odesolve(ODE3, t0..3, Y0)";
feval(symengine, S3)

[Walter Roberson]

%a bit further
Symlist = @(LIST) feval(symengine, 'DOM_LIST@op', LIST);
Symassign = @(VAR, VALUE) feval(symengine, '_assign', VAR, VALUE);

syms y(t)
Dy = diff(y);
D2y = diff(Dy);
IVP3 = Symassign('IVP3', Symlist([D2y == t*Dy, y(0)==0, Dy(0)==2]));
fields3 = Symassign('fields3', Symlist([y, Dy]));

%the [] is important here; it keeps the return value from being the ODE
%and so keeps the MATLAB level from trying to execute the ODE
S3a = "[ODE3, t0, Y0] := [numeric::ode2vectorfield(IVP3, fields3)]; [];";
S3b = "numeric::odesolve(ODE3, t0..3, Y0)";
feval(symengine, S3a);
feval(symengine, S3b)

[Michael Croucher]

This is nice! There's so much useful functionality deep within the MuPAD engine that's not yet fully exposed to MATLAB. Reminds me of an approach I took to vpa QR factorisation several years ago http://walkingrandomly.com/?p=4776 (I think vpa QR factorisation is now fully supported though so none of that is needed)

[Michael Croucher]

Your second code worked for me on Windows/MATLAB 2019a

ans =

[ 70.783509169961080291027822815851, 180.03426260104362710023091349115]

but fails on WindowsMATLAB 2020b

Error using mupadengine/feval_internal
Invalid left side.
Error in mupadengine/feval
Error in symscript (line 13)
feval(symengine, S3a);

[Walter Roberson]

Hmmm, some probing indicates that numeric::odesolve and numeric::ode2vectorfield do not exist in R2020b -- at least not under those names.

I experimented with going back to R2019b and generating a MuPAD notebook containing those functions, and ran the function to convert MuPAD notebooks to LiveScripts -- but unfortunately it was not able to convert those calls.

Hah, I just discovered some earlier work I had posted that did more towards the conversions (but does not solve the R2020b problem); see https://www.mathworks.com/matlabcentral/answers/309163-can-i-use-numeric-odesolve-as-a-replacement-to-ode45#answer_240756

[Walter Roberson]

I have asked Mathworks about the R2020b situation.

[FastCar]

Any luck?

[Walter Roberson]

Mathworks said that those routines are gone and will not be back.

Answer 2

[John D'Errico]

Short answer - no.

Long answer - nnnnoooooo.

ODE45 uses only singles or doubles. Nothing symbolic. vpa is a symbolic tool.

[Michael Croucher]

Thanks John. I thought as much but wanted to check.

[John D'Errico]

I suppose it may happen one day. For example, vpaintegral appeared a few years ago. They could provide a tool to do the same for an ODE. But I don't see it happening soon.

Answer 3

[Steven Lord]

What types of differential equations are you trying to solve that you need to use arbitrary precision arithmetic?

If you're trying to solve the ODE symbolically rather than numerically, take a look at the dsolve function in Symbolic Math Toolbox.

If you have computed a symbolic equation for the right-hand side of your system of ODEs, you can convert that into a form that the numeric ODE solvers in MATLAB (like ode45) can accept that operates on double arrays using odeFunction also in Symbolic Math Toolbox.

This documentation page may be of interest.