Substitute the nonlinearity in the linear system

Hello everyone! I have a question about possibility of substituting some nonlinear function in system of lode’s.

For example, I have a system:

iN=1;%line
jN=2;%row 
Nonlinearity=@(x) x^2;
[T,X]=ode45(@(t,x) odesys(t,x,Nonlinearity,iN,jN), [0 2],[0 0.1]

function RPF=odesys(t,x,Nonlinearity, iN, jN)
  A=[0 1;0 -1];
  B=[0 1];
  A(iN,jN)=A(iN,jN)*Nonlinearity(x(jN))/x(jN);
  RPF=A*x+B;
end

Obviously this works only if x(jN) does’nt equal zero during the integration progress, and I may suppose that this will produce some calculating errors because of dividing.

So I want to know is there better way to do this.

Thank you in advance!

Answer 1

[Walter Roberson]

If you do encounter a division by 0 during calculations, then odesys will end up returning a value that includes either infinite or nan (or both). That will typically not be noticed immediately, but instead the next call to odesys will very likely end up returning mostly or all nan. Everything will fall apart from there.

There are narrow circumstances under which ode45() will notice the big value and recover smoothly, just saying "Oh, I guess that step failed" -- but for that to happen, you have to get lucky enough that the division by 0 was on the very last of the 6 odefun calls per iteration that ode45 makes. Much more typical is that you just start going haywire and figuring out why can be a pain.

So... What you should can is a test such as

RPF(~isfinite(RPF)) = 1e100;   %some big value that is going to be far from normal values

This will eliminate some of the problems, allowing for a cleaner recover.

You should also, however, add an event function through the options structure, and the event function should be marked as Terminal, and Both Directions, and the value to test should be x(2) . An event function such as that will go as close to the border where x(2) becomes 0 as it can, and then will exit ode45 early. You can tell that the exit was early by checking the last time returned.

This particular exception is probably not "removable": if x(jN) becomes 0 because it must due to the initial conditions, then there is not much you can do except quit. There are a lot of other cases where when the exception is detected, that you want to change the boundary conditions. I recommend Mathwork's example ballode which shows how to use events to cause a ball to "bounce" when it hits the floor (it reverses the velocity and then restarts ode45() )

[Ivan Khomich]

Thanks for your answer, Walter!

So, there is no “good” ways to get the nonlinear system from linear, as you’ve mentioned only the ways to avoid the errors?

[Walter Roberson]

I am not familiar with Control Systems theory on "linearizing" a system; I do not know whether that will help.

If you have a pole in your system and that pole is in the area you have to travel through, then you are going to have a problem.

What would you hope would happen if x(jN) becomes 0 due to the natural evolution of the system?

Note:

variable-step solvers such as ode45() work by hypothesizing that a step is valid, and then making some calculations at positions carefully positioned relative to the location (mathematical theory is used about where the positions should be chosen), coming up with a composite result based upon the 6 different locations (in the case of ode45()), and then testing to see whether the change in value is too large relative to current conditions. If it is not too large of an calculated error, then the new proposed position is accepted. If the change is too large, then a shorter step is proposed and it tries again. The step is continually shrunk until either a movement succeeds or the step size gets smaller than the configured minimum, in which case the integration encountered a singularity.

Because of the way that hypothetical positions are calculated for, in the case where the step is too big, then it can happen that a position is evaluated for which the equations are not valid. A point could be tested that cannot actually be reached considering all of the different portions of the ode, but testing at that point can be important for determining which points are valid.

Thus in your situation, you could encounter a case where it proposes x(jN) as 0, perhaps just while testing at the 6 different locations per iteration. The fact that x(jN) becomes 0 does not mean that you always need to immediately quit integration: instead it can work to just return a "big" number to convince ode45 that the step failed and to try again with a shorter step. The event function is for the case where integration must pass through the boundary condition -- the bouncing ball must hit the ground at some point, for example.

But if you have the equivalent of that situation, where x(jN) must become 0 along the way because of the equations, then you want to cleanly exit integration in that case, which you would do using the event function mechanism.

... But in such a situation what would you like ode45() to do? Nonlinear systems can encounter singularities, especially if they have a division by a number that can be 0.

[Ivan Khomich]

Seems, like it is my fault that I phrased my question wrong.

I’ll try it again. Dividing by zeros is the problem that occurs when I create the nonlinear system of differential equations as I’ve shown in my example.

So the real question is: How to create nonlinear system of equations

F=[Nonlin(x(2)); -x(2)+1]

from the matrix of the linear system and anonymous function “Nonlinearity”?

[Walter Roberson]

What you already posted is fine for the case where no singularity can be encountered.

In the case where a singularity can be encountered, what would you like to have happen if the singularity is indeed encountered?

Is the question about having a "black box" function that might have singularities, but at unknown places? If so then the RPF(~isfinite(RPF)) technique avoids accidental singularities during probing of point values, and you can use the nonlinear function inside an ode event function to detect cases where the singularity is inherent.

[Ivan Khomich]

No, the function Nonlin(x) is continuous and twice differentiable, so in the original task there is no singularities. Singularities can appear if the system is created as I show in first message.

And I want to find another way to automatically create this system.

[Walter Roberson]

If your nonlinear function is continuous and twice differentiable, then if you have the Symbolic Toolbox, you could probably construct the forms symbolically and then use odeFunction() to create the anonymous function.

This will not work without a bit of twisting if the function includes min() or max() or mod() : the symbolic versions of those usually do not do what is desired. (On the other hand, those are not differentiable.)

If the nonlinear function has iteration based upon getting down to some particular absolute or relative error, then there could be a problem using symbolic variables with it.

If the nonlinear function compares the inputs to values then a bit of redesign could be needed (but such cases are typically not twice differentiable unless they are carefully constructed.)

[Ivan Khomich]

The functions I need to consider are “good”, so the only problem is constructing the system.

Answer 2

[Ivan Khomich]

I solve this problem by Symbolic Toolbox, as Walter Roberson suggested, so if anyone someday will need this I post my solution:

function fuf
Nonlin=@(x) x^2;
SysOrder=2;
iN=1;
jN=2;
syms xs [1 SysOrder]
A=[0 1;0 -1];
B=[0 ;1];
LinearSystem=A*xs.'+B;
LinearSystem(iN)=subs(LinearSystem(iN), xs(jN),Nonlin(xs(jN)));


[~,X]=ode45(@(t,x) odesys(t,x,LinearSystem, SysOrder), [0 2], [0 0]);

function RPF= odesys(t,x,LinearSystem,SysOrder)
syms xs [1 SysOrder]
RPF=double(subs(LinearSystem, xs(1:SysOrder).', x(1:SysOrder)));

But question partially still open: Is there another way to do this, because of time that code needs to operate with symbolic variables is extremely high(For example, my main programm without this novelty requires 0.5-1s, and 190s with it).

[Walter Roberson]

You should read the documentation for odeFunction and follow the flow in the first example, to come up with a numeric anonymous function to pass to ode45() . Doing symbolic calculations within your ode function is always going to be notably slower than numeric work, and so should be avoided unless there happen to be functions that are only implemented through the symbolic toolbox, or if you need higher precision or higher value range than numeric can give you

[Ivan Khomich]

Thank you very much, Walter!

I tried to do what I‘ve done in last example, but before the ode45 function and get the success in the mobile version. So, I’ll try this later on PC and expect it will work faster.

function fuf
    Nonlin=@(x) x^2;
    SysOrder=2;
    iN=1;
    jN=2;
    syms xs [1 SysOrder]
    A=[0 1;0 -1];
    B=[0 ;1];
    LinearSystem=A*xs.'+B;
    LinearSystem(iN)=subs(LinearSystem(iN), xs(jN),Nonlin(xs(jN)));
    NonlinearSystem=@(x) double(subs(LinearSystem, xs(1:SysOrder).', x(1:SysOrder)));
    [~,X]=ode45(@(t,x) odesys(t,x, NonlinearSystem, SysOrder), [0 2], [0 0]);
    [~,X1]=ode45(@(t,x) [(x(2))^2;1-x(2)],[0 2],[0 0]);
    sum((X-X1).^2)

 function RPF= odesys(t,x, NonlinearSystem,SysOrder)
     RPF= NonlinearSystem(x);

[Ivan Khomich]

The new attempt faster three times, but it still requires a lot of time for calculations(~50 s).

I don't really understand why this happen, beacuse I do not see the difference between numeric system created from symbolic and numeric system handwritten.

If it will not bother you, can you explain why this happen and maybe suggest other way to approach the problem?