Convert the symbolic expression into the expression which can be put on the RHS of an ODE.

3 views (last 30 days)
I have a 2x2 ODE dx1dt = f[1];dx2dt = f[2](see Maple Code at the end for f[1], f[2]). The sensitivity equation is given as dSdt(i,j) = res(i,j) by vector form (corresponding element in matices of equation 5 (time derivative) and 7 in Maple code) for dx3dt~dx8dt appending with state equation for dx1dt and dx2dt, when a=1,b=0,c=1, i.e.,
dx1dt = x2;dx2dt = -sinx1-x2;
dx3dt = x4;dx4dt = -x3cosx1-x2-x4;
dx5dt = x6;dx6dt = -x5cosx1-x6-x2cosx1;
dx7dt = x8; dx8dt = -cos(x1)x7-x8-sinx1; -------------Formula (1)
Matlab code for same function as in Maple code is
nx = 2; % # of state variables
nu = 3; % # of model pars
nx_t = nx + nx*nu;
syms(sym('x', [1, nx_t])) % syms x1 x2 …x8
syms a b c
f(1) = x2
f(2) = -c*sin(x1)-(a+b*cos(x1))*x2
pfpx = jacobian(f,[x1, x2])
pfpw = jacobian(f,[a b c])
S = [x3 x5 x7; x4 x6 x8]
pfpx = subs(pfpx,[a,b,c],[1,0,1])
res = pfpx*S + pfpw
I want to solve the sensitivity equation in Matlab, I could type the function by Formula(1) on the RHS manually as,
function dxdt = dyneqn(t,x)
dxdt = zeros(8,1)
dxdt(1) = x(2);
dxdt(2) = -sin(x(1))-x(2);
dxdt(3) = x(4);
dxdt(4) = -x(3)*cos(x(1))-x(2)-x(4);
dxdt(5) = x(6);
dxdt(6) = -x(5)*cos(x(1))-x(6)-x(2)*cos(x(1));
dxdt(7) = x(8);
dxdt(8) = -cos(x(1))*x(7)-x(8)-sin(x(1));
end
If allowing user to input different f(1) and f(2) expressions, is there an efficient way to auto-generate ODE function ‘dyneqn’ as shown above based on ‘res’ (= pfpx*S + pfpw)? The difficult part is the state variables in ‘res’ are symbolic variables x1,…,x8, in order to express it in ‘dyneqn’, needs to be like x(1), x(2),…,x(8) on RHS of ODE, and one can not define symbolic variables like x(1),…,x(8) in Matlab. This problem is tricky and really appreciate in advance!

Answers (1)

Hainan Wang
Hainan Wang on 12 Jun 2021
I work this out with odefunction. Guess I will answer this myself.
clear all
clc
%% Part1.Define symbolic variables
nx = 2; % # of state variables
nu = 3; % # of model pars
nx_t = nx + nx*nu;
% 1.1 Define Variables as symfun
variableName = sprintfc('x%d(t)', 1:nx_t);
syms(variableName{:})
vars = str2sym(sprintfc('x%d(t)', 1:nx_t));
syms a b c
% vars =[x1(t),x2(t)] 1x2 sym
% x1 % 1x1 symfun
%1.2 Define Variables as sym
% syms a b c
% syms(sym('x%d', [1, nx_t]))
% vars = sym('x',[1 nx])
% % x1 % 1x1 sym
% vars = [x1, x2] 1x2 sym
%% Part2.Build sensitivity function
f(1,1) = x2;
f(2,1) = -c*sin(x1)-(a+b*cos(x1))*x2;
pfpx = Jacob_fun(f,nx);
% pfpx = jacobian(f,sym('x',[1, nx]))
pfpw = jacobian(f,[a b c]);
S = formula([x3 x5 x7; x4 x6 x8]);
pfpx = subs(pfpx,[a,b,c],[1,0,1]);
res = pfpx*S + pfpw;
f(2) = subs(f(2),[a,b,c],[1,0,1]);
for i = nx+1: nx_t % Be advised, size f change
f(i,1) = res(i-nx);
end
M = sym(diag(ones(1,4)));
odefun = odeFunction(f,vars);
initiConditions = [1 1 zeros(1,nx*nu)]';
%% Part3.Visualization
[t,z] = ode15s(odefun,[0 10],initiConditions)
subplot(2,1,1)
plot(t,z(:,[3 5 7]))
subplot(2,1,2)
plot(t,z(:,[4 6 8]))

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!