Error using odearguments with a symbolic system of differential equations

Hello Everyone,

I got the following errors:

Error using odearguments (line 113)
Inputs must be floats, namely single or double.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in SxMatrix_grinding (line 11)
[t1, sol1] = ode45('Sxd1', tspan1, IC1);

My function is:

function output=Sxd1(t, X)
syms u x1 x2 k1 k2 k3 k4 k5 k1_ k2_ k3_ k4_ k5_ x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14

u=0.0254;

var=[x1;x2];
k=[k1_;k2_;k3_;k4_; k5_];

f=[k1*k1_*x2;
    (u-x2)/(k2*k2_+k3*k3_*x1)];

dfdk=jacobian(f,k);
dfdx=jacobian(f,var);

Sx = [x3 x5 x7 x9 x11; x4 x6 x8 x10 x12];

Sxd=dfdx*Sx+dfdk;

x1=X(1);
x2=X(2);
x3=X(3);
x4=X(4);
x5=X(5);
x6=X(6);
x7=X(7);
x8=X(8);
x9=X(9);
x10=X(10);
x11=X(11);
x12=X(12);

output=[k1*k1_*x2;
    (u-x2)/(k2*k2_+k3*k3_*x1);
    Sxd(:)];
output=subs(output, {k1, k2, k3, k4, k1_, k2_, k3_, k4_}, {220, 0.7463, 0.00216, 88000, 1, 1, 1, 1});

end

and my main code is:

tspan1 =[0 9.5];

%options = odeset('RelTol', 1e-36, 'AbsTol', 1e-50);
IC1=zeros(1,12); IC1(11)=1;
[t1, sol1] = ode45('Sxd1', tspan1, IC1);

I am having issues understaing the reason of the error I got. Any help is really appreciated!

Thank you in advance!

[Walter Roberson]

dfdk=jacobian(f,k);

It is a waste of time to do that calculation every time. Do it once before-hand, and use odeFunction() or matlabFunction() to convert to a function handle that works on numeric values.

[Gabriele Galli]

Hello Walter,

thank you so much for your reply.

I read about odeFunction() and matlabFunction().

I read the documentations of both the functions but I cannot fully understand the differences between them.

Moreover, I try to use odeFunction in the following way:

vars=[x1 x2 x3 x5 x7 x9 x11 x4 x6 x8 x10 x12]
output = odeFunction(output, vars)

but I got the following errors:

Error using mupadengine/feval_internal
Variables must be specified by univariate function calls.

Error in sym/odeFunction (line 118)
  A = feval_internal(symengine, 'daetools::odeFunction', expr, vars, params{:});

Error in untitled4 (line 26)
output = odeFunction(output, vars)

I am so sorry but I am not familiar at all with symbolic tool.

Thank you so much.

[Walter Roberson]

syms u x1 x2 k1 k2 k3 k4 k5 k1_ k2_ k3_ k4_ k5_ x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14

at that point, x1 through x14 are symbolic

var=[x1;x2];
f=[k1*k1_*x2;
    (u-x2)/(k2*k2_+k3*k3_*x1)];
...
Sx = [x3 x5 x7 x9 x11; x4 x6 x8 x10 x12];

Those use the symbolic x1 through x12 values.

Sxd=dfdx*Sx+dfdk;

Still symbolic.

x1=X(1);
x2=X(2);
etc

Okay now, at the MATLAB level, x1 through x12 are assigned numeric values. But only at the MATLAB level. This does not affect Sxd, which still uses the symbolic variables.

The situation is exactly like

A = 1
B = A + 1
A = 2

B does not suddenly become 3. B is not a formula: at the time that B was assigned to, the current value of A was copied and used in the computation to create B, and B then promptly forgets that it knew anything about the source variable A.

Just the same way if you have

syms A
B = A + 1
A = 2

then B does not suddenly become 3. At the time that B was assigned to, the current value of A was copied and used in the computation to create B. That current value is sym('A'), a reference to a symbol that lives inside the symbolic engine. B promptly forgets that it knew anything about where that sym('A') came from. B is not "reference to variable A, dereference, add 1": B is sym('A')+1 and no change to MATLAB variable A will affect B unless you tell MATLAB to force it to.

output=[k1*k1_*x2;
    (u-x2)/(k2*k2_+k3*k3_*x1);
    Sxd(:)];

At that point, the x1 and x2 that you can see in the first two lines of output have become the input numeric variables. However, Sxd has not been affected, and still contains vectors of symbolic variables including sym('x1') and sym('x2')

output=subs(output, {k1, k2, k3, k4, k1_, k2_, k3_, k4_}, {220, 0.7463, 0.00216, 88000, 1, 1, 1, 1});

Okay, you have explicitly told MATLAB to look inside output to find references to the variables indicated by those k* and k*_ and replace them with specific numeric values. That is generally a good thing to do.

But you are not telling it to replace any of the symbolic x* values, so output still potentially contains references to symbolic sym('x1') through sym('x14') . And that is why the double() fails.

[Gabriele Galli]

Thank you so much for the clarification!

Now, I understand how it works and the reason of the error.

The last point I did not get is how can I replace the symbolic values (x1 through x14) inside my output value.

Can you help me with this?

Thank you!

[Walter Roberson]

output = subs(output);

[Gabriele Galli]

Thank you,

I tried to add to the function

...
    output=[k1*k1_*x2;
    (u-x2)/(k2*k2_+k3*k3_*x1);
    Sxd(:)];
output=subs(output, {k1, k2, k3, k4, k1_, k2_, k3_, k4_}, {220, 0.7463, 0.00216, 88000, 1, 1, 1, 1});

output=subs(output);
end

but unfortunately I am still getting the same error:

Error using odearguments (line 113)
Inputs must be floats, namely single or double.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in SxMatrix_grinding (line 11)
[t1, sol1] = ode45('Sxd1', tspan1, IC1);

[Walter Roberson]

Try with @Sxd1 instead of 'Sxd1'

Double check class() tspan1 and IC1

[Gabriele Galli]

I am still getting the same error. This is the way I called the function:

tspan1 =[0 9.5];
IC1=zeros(1,12); IC1(11)=1;
[t1, sol1] = ode45(@(t,X) Sxd1(t,X), tspan1, IC1);

Plus, tspan and IC1 looks right to me.

Answer 1

[Steven Lord]

Your ODE function can use symbolic calculations internally but it must return a double or single array to the ODE solver. Call double on your output immediately before the end of your ODE function.

Alternately if you want to solve it symbolically (and if you want to use such stringent tolerances you may) use the functions included in Symbolic Math Toolbox for solving differential equations like dsolve.

[Walter Roberson]

Note that even if your ode function gets down to a symbolic number line sym('0.342379235323259883251') then that is not good enough for ode45(): the output of the function needs to be double() or single()

[Gabriele Galli]

Hello Steven,

thank you for the help, I have already tried using

output=double(output)

but still doesn't work, I got the following error:

Error using symengine
Unable to convert expression into double array.

Error in sym/double (line 698)
        Xstr = mupadmex('symobj::double', S.s, 0);

Error in Sxd1 (line 37)
output=double(output)

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in SxMatrix_grinding (line 11)
[t1, sol1] = ode45('Sxd1', tspan1, IC1);

Concerning the tolerance, it is ok using a less stringent tolerance as well.

Thank you!

[Steven Lord]

Set an error breakpoint and run your code. When MATLAB stops, select the workspace of the Sxd1 function and look at the contents of the variable output. Does it contain any symbolic variables? [It does.] If it does, substitute values for all those symbolic variables into output then try converting it to double again.

[Walter Roberson]

We would need to see your current code to say more.

But I repeat my advice from https://www.mathworks.com/matlabcentral/answers/596998-error-using-odearguments-with-a-symbolic-system-of-differential-equations#comment_1016971 : Do not do that symbolic calculation each time. Do the symbolic calculation once and use odeFunction or matlabFunction to generate a function handle for you.