Calculating the Jacobian -- How to Convert an Anonymous Function to a Symbolic Function

I am currently working on a complicated system of ~30 ODE's. To solve the system I have to use a stiff solver and calculate the jacobian of the system. I have written a successful series of functions that accomplishes this, however it requires that I define the ODE system in the "jacobian-calculating" function. Because of the complexity of my system I would prefer to create a stand-alone "model-function" that can be called by either my ode solver or my jacobian-calculating function.

I have included minimal examples of my code below. I believe my problem is with the line:

-----------------------------------------------------

%Convert the "Numeric" System to a Symbolic System

SymbolicSystem = sym(NumericSystem);

-----------------------------------------------------

where "NumericSystem" is the handle for the anonymous function of my ode system. This line is contained in my function "testjacobian".

Does anyone know how to convert an anonymous function into a symbolic function? Or is there a better strategy for what I'm trying to do?

I do know that my code seems a bit contrived for this simple system, but this is because I have tried to capture the essential elements of the code used to solve my more complex system.

Any help would be greatly appreciated, Thanks in advance,

Laura

-----------------------------------------------------------------------------------

CODE:

-----------------------------------------------------------------------------------

function testdriver

%Parameters

a = 1;

b = 2;

c = 3;

d = 4;

param = [a b c d];

%Initial Conditions

start1 = 0;

start2 = 1;

init = [start1, start2];

%Calculate the solution

[t,x] = testsolver(param,init);

%Dynamics Plot

%Call a function that plots the data

% .... %

end

-----------------------------------------------------------------------------------

function [t,x]=testsolver(param,initial)

%Integration Domain

tspan = [0,10];

%Initial Conditions

x1init = initial(1);

x2init = initial(2);

xinit = [x1init,x2init];

%Construct the Scalar-Valued, Numeric, Jacobian Function

J = testjacobian(param);

%Convert the function to a Vector-Valued, Numeric Function

CalcJac = @(t,x) call_on(J,x);

%Define Integration Options -- Use the Jacobian

opts1 = odeset('Jacobian',CalcJac); 

%Integrate the System

[t,x] = ode15s(@(t,x)testmodel(t,x,param),tspan,xinit,opts1);

end

function r = call_on(func, array)

%This Function Converts a Scalar-Valued Function to a Vector-Valued Function

      vector = num2cell(array);

      r = func(vector{:});

end

-----------------------------------------------------------------------------------

function xp = testmodel(t,x,param)

%Define Parameters

A = param(1);

B = param(2);

C = param(3);

D = param(4);

%Assign x values to Variable Names

thing1 = x(1);

thing2 = x(2);

%Initialize the Derivatives to Zero

xp = zeros(2,1);

%System of ODE's

xp(1) = A*thing1^2 + B*thing2;

xp(2) = C*thing1 + D;

end

-----------------------------------------------------------------------------------

function Jacob=testjacobian(param)

%Initialize Symbolic Variables

syms tsym y1 y2

y = [y1 y2];

%Using "param" Construct a (t,x)-Dependent System

NumericSystem = @(t,x)testmodel(t,x,param);

%Convert the "Numeric" System to a Symbolic System

SymbolicSystem = sym(NumericSystem);

%Calculate the Jacobian of the Symbolic System

SymbolicJacobian = jacobian(SymbolicSystem',y);

%Convert the Symbolic Function to a "Numeric" Function Dependent on y

Jacob = matlabFunction(SymbolicJacobian,'var',y);

end

[Laura]

Edited Original Question:

Original Text contained the line: " [t,x,TE,YE,IE] = testsolver(param,init);"

Because the minimal code provided does not contain any event definitions the provided code has been updates to read: " [t,x] = testsolver(param,init);"

Answer 1

[Walter Roberson]

Replace

NumericSystem = @(t,x)testmodel(t,x,param);
%Convert the "Numeric" System to a Symbolic System
SymbolicSystem = sym(NumericSystem);

with

syms t
SymbolicSystem = testmodel(t, y, param);

You are thinking of testmodel as numeric, but really at this point you are dealing in formula to build the symbolic jacobian, so it is fine to pass symbolic variables to it. You would have difficulty if the actual code in testmodel used "if" or comparisons based upon any parameter that you pass in symbolically, but such systems could be written in terms of Heaviside (but watch out for what Heaviside(0) means.)

Remember to simplify() your symbolic jacobian before turning it into an anonymous function, just for efficiency.

[Laura]

Thank you Walter, I will follow up once I've had a chance to try this out.

[Laura]

Walter, thank you again for taking the time to answer my question. Your solution worked with one modification. My goal was to structure testmodel in such a way that it could be called (in testdriver) by both testjacobian and testsolver.

ode15s requires (as I understand it) that the derivative vector(xp) in the model be initialized to the zero vector.

Initializing the derivative vector (xp) to zeros in testmodel causes it to become type "double" which causes problems when trying to calculate the symbolic Jacobian.

Using your suggestion above and changing

>xp = zeros(2,1);

to

>if isa(t,'double')
xp = zeros(2,1);
end

allowed the codes to run the way I'd hoped.

Thank you for your advice,

Laura