How can I plug the symbol vector into function.m file for ode solving?

I want to solve the problem where a stage is that I need to use ode23 or ode45 to solve the time-varying variables. For instance, the varibales is nX1 space vector. So I need to build the nX1 differential equations like that

function f = ODE(t,y,Q,R,r)

f =[2*y(1) - 10*exp(-2*t) - 22*y(3) - (y(1)*(y(1)/2 + 25*y(2) - (51*y(3))/2))/(2*(t - 15)) - (25*y(2)*(y(1)/2 + 25*y(2) - (51*y(3))/2))/(t - 15) + (51*y(3)*(y(1)/2 + 25*y(2) - (51*y(3))/2))/(2*(t - 15))

51*y(2) - 60*y(3) - 11*y(6) - (y(1)*(y(2)/2 + 25*y(5) - (51*y(6))/2))/(2*(t - 15)) - (25*y(2)*(y(2)/2 + 25*y(5) - (51*y(6))/2))/(t - 15) + (51*y(3)*(y(2)/2 + 25*y(5) - (51*y(6))/2))/(2*(t - 15))

11*y(3) - 10*y(4) - 11*y(8) - (y(1)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/(2*(t - 15)) - (25*y(2)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/(t - 15) + (51*y(3)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/(2*(t - 15))...]

Traditionally, I can solve it by using [t,y]=ode45(@ODE,tspan,InitialCondition) to find the curve of all of the variables.

My point is: Can I let the differnential equations be automatically solved as soon as I get the symbolic expression of matrix f. I have think about it for serval days already.

Thanks.

[Stephan]

Can you give an example on how the input looks like and what you expect to get out? Is it always the same system of ode's with just other coefficients?

[Zhichao Shen]

For example, the variables are y1, y2, ....y14. And the input include the state matrices like A, B, Q(t), R(t), r(t). Then I already get the differnential equations about the 14 variables. Like what I said before, the state matrices Q,R,r is dependent on time, so if I change the expression of one of them, I will derive the differnent ODEs.

So I need to change the expression of ODEs everytime I change one of them(Q(t),R(t),r(t)). The problem is located in control theory. The detail brach is functional approach for solving linear tracking problem. The output is the result of y(t) by ode23 or ode45, which is discrete on time since it's totally impossible to find the close-form expression of these variables.

Not like what you said, the coefficients is time-varying.

This is the ODEs of one of the examples I can show. It's a 14 variables problem. I can get the content in f by given A,B,Q,R,r as pretreatment of a code file I finished. The expression is a little comlicated. So I don't hope I need to copy and paste every time.

function f = ODE(t,y,Q,R,r)

f =[2*y(1) - 22*y(3) - Q(1, 1) + (y(1)*(y(1)/2 + 25*y(2) - (51*y(3))/2))/(2*R) + (25*y(2)*(y(1)/2 + 25*y(2) - (51*y(3))/2))/R - (51*y(3)*(y(1)/2 + 25*y(2) - (51*y(3))/2))/(2*R)

51*y(2) - 60*y(3) - 11*y(6) - Q(2, 1) + (y(1)*(y(2)/2 + 25*y(5) - (51*y(6))/2))/(2*R) + (25*y(2)*(y(2)/2 + 25*y(5) - (51*y(6))/2))/R - (51*y(3)*(y(2)/2 + 25*y(5) - (51*y(6))/2))/(2*R)

11*y(3) - 10*y(4) - 11*y(8) - Q(3, 1) + (y(1)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/(2*R) + (25*y(2)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/R - (51*y(3)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/(2*R)

11*y(4) - 11*y(9) - Q(4, 1) + (y(1)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/(2*R) + (25*y(2)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/R - (51*y(3)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/(2*R)

100*y(5) - 120*y(6) - Q(2, 2) + (y(2)*(y(2)/2 + 25*y(5) - (51*y(6))/2))/(2*R) + (25*y(5)*(y(2)/2 + 25*y(5) - (51*y(6))/2))/R - (51*y(6)*(y(2)/2 + 25*y(5) - (51*y(6))/2))/(2*R)

60*y(6) - 10*y(7) - 60*y(8) - Q(3, 2) + (y(2)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/(2*R) + (25*y(5)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/R - (51*y(6)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/(2*R)

60*y(7) - 60*y(9) - Q(4, 2) + (y(2)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/(2*R) + (25*y(5)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/R - (51*y(6)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/(2*R)

20*y(8) - 20*y(9) - Q(3, 3) + (y(3)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/(2*R) + (25*y(6)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/R - (51*y(8)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/(2*R)

20*y(9) - 10*y(10) - Q(4, 3) + (y(3)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/(2*R) + (25*y(6)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/R - (51*y(8)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/(2*R)

20*y(10) - Q(4, 4) + (y(4)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/(2*R) + (25*y(7)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/R - (51*y(9)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/(2*R)

r(1)*Q(1, 1) + r(2)*Q(1, 2) + r(3)*Q(1, 3) + r(4)*Q(1, 4) + y(11)*((y(1)/2 + 25*y(2) - (51*y(3))/2)/(2*R) + 1) - y(13)*((51*(y(1)/2 + 25*y(2) - (51*y(3))/2))/(2*R) + 11) + (25*y(12)*(y(1)/2 + 25*y(2) - (51*y(3))/2))/R

r(1)*Q(2, 1) + r(2)*Q(2, 2) + r(3)*Q(2, 3) + r(4)*Q(2, 4) + y(12)*((25*(y(2)/2 + 25*y(5) - (51*y(6))/2))/R + 50) - y(13)*((51*(y(2)/2 + 25*y(5) - (51*y(6))/2))/(2*R) + 60) + (y(11)*(y(2)/2 + 25*y(5) - (51*y(6))/2))/(2*R)

r(1)*Q(3, 1) - 10*y(14) + r(2)*Q(3, 2) + r(3)*Q(3, 3) + r(4)*Q(3, 4) - y(13)*((51*(y(3)/2 + 25*y(6) - (51*y(8))/2))/(2*R) - 10) + (y(11)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/(2*R) + (25*y(12)*(y(3)/2 + 25*y(6) - (51*y(8))/2))/R

10*y(14) + r(1)*Q(4, 1) + r(2)*Q(4, 2) + r(3)*Q(4, 3) + r(4)*Q(4, 4) + (y(11)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/(2*R) + (25*y(12)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/R - (51*y(13)*(y(4)/2 + 25*y(7) - (51*y(9))/2))/(2*R)];

end

[Stephan]

If the shown equations already are the result of your coding, then the way stated in Ameer's answer is the correct one. Additionally you will need

matlabFunction

to get the result of

odeToVectorField

suitable for

ode45

.

[Walter Roberson]

Not matlabFunction: use the workflow for odeFunction()

[Stephan]

Hi Walter, whats the advantage of odeFunction? I only know it from usage with algebraic systems from the documentation.

In the doc for odeToVectorField they use matlabFunction:

https://de.mathworks.com/help/symbolic/odetovectorfield.html#bvlvtw5-1

Can you tell more?

[Zhichao Shen]

I have same question. I have saw the doc of odeFunction. It seems that this code is situable for variables, which form like x(t) y(t) and so on.

[Walter Roberson]

When you build up ode equations symbolically you use symbolic functions x(t) and so on. The workflow converts the function and derivatives references into variable references so that the boundary condition vector can be indexed, and odeFunction knows to package in the form of f(t,x, additionals)

[Stephan]

Thank you.

Answer 1

[Ameer Hamza]

Are you manually pasting your symboling expression in the ODE() function? You can use odetovectorfield(): https://www.mathworks.com/help/symbolic/odetovectorfield.html to solve the equations with ode45(). See the example "Solve Higher-Order Differential Equation Numerically" on the documentation page of odetovectorfield().

[Zhichao Shen]

Yes I was pasted the content in ODE manually. So I think it's unconvinent if I need to more steps which maybe based on the differnent inputs in future. I will try the code you suggested. Thanks.

[Zhichao Shen]

Hi, Ameer

I already let my code be ccompatible with this method. I have one more question about this question, how can I set

syms y1(t) y2(t) ... yn(t)

I have setted it manually.

But I hope I can solve any problem for any given n.

Thanks.

[Ameer Hamza]

Yes, it is possible to create symbolic functions like this automatically. Read some detail about this in comments here: https://www.mathworks.com/matlabcentral/answers/652078-how-to-define-recursively-symbolic-function-for-differential-equation-system. In short, you can use something like this

n = 10;
syms(compose('y%d(t)', 1:n))

[Zhichao Shen]

Thanks a lot!!! You are so nice！

[Ameer Hamza]

I am glad to be of help! :)

[Zhichao Shen]

Hi, Ameer

I still don't know how to create a vector y=[y1(t) y2(t) ... yn(t)] after I build these elements afters yms(compose('y%d(t)', 1:n))

I have looked the link you sent.

[Ameer Hamza]

Walter's comment mentions about this

syms t
n = 10;
y = cellfun(@(S) symfun(S, t), compose('y%d(t)', 1:n), 'uniform', 0)
A = [A{:}]

[Zhichao Shen]

Hi, Ameer

I was didn't find a mistake under this approach.

The mistake is that the sequence of first one and the second one differential equation are switched. I don't know what happed excatly result the problem.

The code and result is like that:

>> s=[diff(y1,t)==2*y2(t)^2;diff(y2,t)==y1(t)^2+y2(t)]

s(t) =

diff(y1(t), t) == 2*y2(t)^2

diff(y2(t), t) == y1(t)^2 + y2(t)

>> odeToVectorField(s)

ans =

Y[2]^2 + Y[1]

2*Y[1]^2

but I want it should be like(I type manually)

2*Y[2]^2

Y[1]^2 + Y[2]

Maybe I can't use odeToVectorField. The form is differnent a little bit.

[Ameer Hamza]

If all of them are first order ODE, then you can get the second output of odeToVectorField and sort it to get the order you want.

syms y1(t) y2(t)

s = [diff(y1,t)==2*y2(t)^2;diff(y2,t)==y1(t)^2+y2(t)];

[V, S] = odeToVectorField(s);
[~, idx] = sort(S);
V = V(idx)

[Walter Roberson]

When you [] expressions that involve symbolic functions, then the result is a symbolic function that returns a vector

syms x(t) y(t)
z = [x y]

z will be a symfun that returns a vector, not a vector of symbolic functions. It is not possible to be able to use () indexing to return a symfun: subsref with () indexing for symfun is defined as function invocation.

You can [] together symfun and symbolic expressions and the resulting symfun will have the same calling sequence as the original symfun

You cannot [] two symfun with different variables or different orders of variables.

My code created cellfun because of this behavior of [] creating a single symfun that returns a vector.

[Ameer Hamza]

Ah! Yes, that's correct. Not sure how OP is going to use it, but the correct is to create a cell array of symbolic functions.

[Zhichao Shen]

Ameer,

The output of your code is also

2*Y[1]^2

Y[2]^2 + Y[1]

You have let the sequence be changed. But the sequence of the elements y1, y2 no change.

[Ameer Hamza]

I don't think there is a good way to change that. But why does this order matter? The solution will still be the same.

Answer 2

[Zhichao Shen]

Ameer and Walter

Until now, I can say I make a breakthrough of the problem by using odeFunction(dydt,vars) where dydt is a series of ODEs that I derived like I showed originally, then vars is the vector of [y1(t), y2(t), ... yn(t)].

Please keep attention that dydt is set of right hand of ODEs, for instance,

dydt=[y2(t);y1(t)*y2(t)]

,which means y1'(t)=y2(t) and y2'(t)=y1(t)*y2(t) as default setting.

Then we can use ode function to solve the problem for any given inputs and n(variables).

All in all, the main and significant contribution is from Ameer Hamza and Walter Roberson. And I am glad Stephan response to me quickly.

If someone have similar question in the future, please see this summary answer.

Thank you, guys!