Convert the symbolic expression into the expression which can be put on the RHS of an ODE.
3 次查看(过去 30 天)
显示 更早的评论
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!
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/615115/image.jpeg)
2 个评论
Walter Roberson
2021-5-11
See matlabFunction() and odeFunction() . See in particular the workflow in the first example to odeFunction()
回答(1 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!