Computing the jacobian of an anonymous function - MATLAB

I'm trying to solve a system of non linear odes in Matlab as follows.

editparams %parameters

Tend = 50;

Nt = 100;

% Define RHS functions

RHS = @(t,x) ModelRHS(t,x,param); %anonymous function defining the set of equations

%Execution-----------------------------------------------------------------

x0 =[0.04;0.75;0.85]; %Initial condition

t = linspace(0,Tend,Nt); %TSPAN

[t x] = ode45(RHS, t, x0);

Now, I need to find the steady state of the system and I'm trying to create a function for this. I thought I'd calculate the steady state using the Jacobian. My equations are in an anonymous function which is defined as f in the code below. However, I realised that jacobian does not work for anonymous functions (or may be there is a way to do with anonymous functions) . so I thought I would convert the anonymous function to a symbolic function and try it. But i still have a difficulty in getting that done. So any help is really appreciated!

function SS = SteadyState(param, E)

f = @(t,x)ModelRHS(t,x,param,E); %set of odes

SymbolicSystem = sym(f); %trying to make the anonymous function symbolic

SymbolicJacobian = jacobian(SymbolicSystem',x); %jacobian 
Jacob = matlabFunction(SymbolicJacobian,x);

end

Also if there is any other way apart from finding the Jacobian, kindly let me know about hat too.

[UserCJ]

This is exactly the same question that I posted here and here.

[Torsten]

So somebody stole your question or is your pseudonym "Ali Raza" ?

[UserCJ]

I don't use a pseudonym btw!

Answer 1

[Torsten]

Why not simply using "fsolve" to calculate steady state ?

For "fsolve" to work, you only have to supply the function ModelRHS - computation of the Jacobian will be done by "fsolve" internally.

[UserCJ]

Thanks, Torsten for your comment. I did like this;

f = @(t,x)ModelRHS(t,x,param);
x0 =[0.04;0.75;0.85];
options = optimoptions('fsolve','Display','iter'); % Option to display output
SS = fsolve(f,x0,options); % Call solver

But returned an error

Not enough input arguments.

Error in @(t,x)ModelRHS(t,x,param)

Error in fsolve (line 242)
fuser = feval(funfcn{3},x,varargin{:});

Caused by:
Failure in initial objective function evaluation. FSOLVE cannot continue.

[Torsten]

f = @(x)ModelRHS(t,x,param);

instead of

f = @(t,x)ModelRHS(t,x,param);

I assume f does not explicitly depend on t (otherwise there would not be a steady-state).

So set t = 0, e.g., before defining f as above.

I assume further that the necessary parameters "param" were set before f is defined.

[UserCJ]

I did the following

t=5
f = @(x)ModelRHS(t,x,param);
x0 =[0.04;0.75;0.85]';

options = optimoptions('fsolve','Display','iter'); % Option to display output
SS = fsolve(f,x0,options); % Call solver

and it worked; and this is what I got.

                                      Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          4         3.61108                          3.64               1
     1          8         0.85106              1           1.55               1
     2          9         0.85106              1           1.55               1
     3         13         0.45359           0.25          0.666            0.25
     4         17       0.0467741          0.625          0.468           0.625
     5         21      0.00116572       0.225104         0.0967            1.56
     6         25     2.55005e-08      0.0238037       0.000606            1.56
     7         29     1.41202e-18    4.15513e-05       4.51e-09            1.56

But I have no idea what these things are, and I'm specifically looking at negative and positive steady state values. Is there anyway that I can get it coded?

[Torsten]

SS contains the steady states for the differential equations you implemented.

So if you type

SS

after the line

SS = fsolve(f,x0,options); % Call solver

you get the values of the solution variables of your differential equations in steady state.

What do you mean by

Is there anyway that I can get it coded?

[UserCJ]

Oh yes, I got it! Thanks Torsen!

What if I define the initial condition using the steady state function, for which the steady state at t=0? Like the code below, I'm defining my initial guess to be the steady state value at t=0.

x0 = SteadyState(param, 0)

Then, I cannot pass x0 to fsolve to get the SS.

editparams %parameters
Tend = 50;
Nt = 100;
% Define RHS functions
RHS = @(t,x) ModelRHS(t,x,param,E); %anonymous function defining the set of equations
%Execution-----------------------------------------------------------------
x0 = SteadyState(param, 0); %Initial condition as the SS at t=0
t = linspace(0,Tend,Nt); %TSPAN
[t x] = ode45(RHS, t, x0);

function SS = SteadyState(param, E)

t=5
f = @(x)ModelRHS(t,x,param,E);
x0 =[0.04;0.75;0.85]';

%options = optimoptions('fsolve','Display','iter'); % Option to display output
SS = fsolve(f,x0); % Call solver
end

[Torsten]

Then, I cannot pass x0 to fsolve to get the SS.

I don't understand what you mean. The x0 is the SS from fsolve if you call

x0 = SteadyState(param, 0)

Please describe again what your problem is.

[UserCJ]

If you notice the SS code chunck, we need to specify x0 to use fsolve to get the steady state.

t=5
f = @(x)ModelRHS(t,x,param);
x0 =[0.04;0.75;0.85]'; %initial condition
%options = optimoptions('fsolve','Display','iter'); % Option to display output
SS = fsolve(f,x0);

But at some point if we want to define the initial condition through the steady state function as given below, how can we do that?

x0 = SteadyState(param, 0);

The reason for asking this is, when we call fsolve, we need to give x0 as

SS = fsolve(f,x0)

How can we give x0 without knowing that?

[Torsten]

If you have knowledge about approximate values for SS, you can use this information to prescribe a value for x0 that is in the neighbourhood of SS. But in principle, x0 can be chosen arbitrarily. fsolve should converge to a solution of RHS=0 independent of x0.

[UserCJ]

Yes, in this steady state calculation and model simulation, I have chosen x0 arbitrarily. However, in some stage I need my model to take it's x0 through the steady state. So I need to call steady state function at t=0 like;

x0 = SteadyState(param, 0);

But the steady state function is created using fsolve and it needs x0 as an input

fsolve(f,x0)

So without chosing x0 arbitrarily, how can we feed x0 to fsolve or is there anyway that I can do these two things?

I'm sorry, I'm new to matlab and still learning! Appreciate your help Torsen!

[Torsten]

But the steady state function is created using fsolve and it needs x0 as an input

No, it does not need your steady-state x0 as input. It needs an arbitrary vector of length 3 as input.

The syntax is the following:

% Calculate steady-state of your model
xfsolve_start = rand(3,1);                  % create random vector for fsolve to start with
x0 = SteadyState(xfsolve_start,param,E);    % calculate steady-state of the model equations
% Prepare time-dependent calculation
Tstart = 0.0;
Tend = 50;
Nt = 100;
% Define RHS functions
RHS = @(t,x) ModelRHS(t,x,param,E);          %anonymous function defining the set of equations
%Execution-----------------------------------------------------------------
t = linspace(Tstart,Tend,Nt);                %TSPAN
[t ,x] = ode45(RHS, t, x0);                  % Start your time-dependent calculation in steady-state

function x0 = SteadyState(xfsolve_start,param,E)
  t = 0;                                       % Chosen arbitrarily ; doesn't influence steady-state solution since RHS does not explicitly depend on t (I hope)
  RHS = @(x)ModelRHS(t,x,param,E);
  x0 = fsolve(RHS,xfsolve_start);
end

[UserCJ]

Thank you very much Torsten! Your comments helped me to understand the whole situation.

Answer 2

[Matt J]

On the File Exchange, there are numerical methods for approximating the Jacobian of a non-symbolic function using finite differences , e.g.,

https://www.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation

You should probably add a post-processing step that smooths the output of ModelRHS() with a smoothing spline, however. ODE functions seem to have a tendency to produce noisy, non-differentiable output, see,

https://www.mathworks.com/help/optim/ug/optimizing-a-simulation-or-ordinary-differential-equation.html#btfb69a