Problem related with "odeToVectorField"

I used this code to create symbolic equation.

a_sym = sym('a', [N, 1], 'real');
a_dot_sym = sym('a_dot', [N, 1], 'real');
dJ_da_sym = sym('dJ_da', [N, 1], 'real');
for j = 1:N
    a_dot_sym(j) = dJ_da_sym(j)*(dJ_da_sym(j) - lambda);
end

I got these equations 3 equation in ode.

-(100000*(4*a1 + 1/25)*(a1 - 99/100)^2 - 292000*a1*(a1 - 1)^2 + 292*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 - 800*a2*((53*a1)/500 - 499947/50000)^3*(a2 - 1)^2 - 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2 + 1200*a3*((53*a2*((53*a1)/5000 - 499947/500000))/500 + 10)^3*(a3 - 1)^2)*(292000*a1*(a1 - 1)^2 - 100000*(4*a1 + 1/25)*(a1 - 99/100)^2 - 292*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 + 800*a2*((53*a1)/500 - 499947/50000)^3*(a2 - 1)^2 + 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2 - 1200*a3*((53*a2*((53*a1)/5000 - 499947/500000))/500 + 10)^3*(a3 - 1)^2 + 1/100)
                                                                   (108000*a1*(a1 - 1)^2 - 100*(8*a2 + 2/25)*((53*a1)/500 - 10)^3*(a2 - 99/100)^2 + 1200*a3*(a3 - 1)^2*((53*(2*a2 + 1/50)*((53*a1)/5000 - 1))/1000 + 10)^3 + 292*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 - 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2)*(108000*a1*(a1 - 1)^2 - 100*(8*a2 + 2/25)*((53*a1)/500 - 10)^3*(a2 - 99/100)^2 + 1200*a3*(a3 - 1)^2*((53*(2*a2 + 1/50)*((53*a1)/5000 - 1))/1000 + 10)^3 + 292*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 - 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2 - 1/100)
                                                                                                                                                                            -(108000*a1*(a1 - 1)^2 + 100*(12*a3 + 3/25)*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 99/100)^2 - 508*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 - 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2)*(508*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 - 100*(12*a3 + 3/25)*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 99/100)^2 - 108000*a1*(a1 - 1)^2 + 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2 + 1/100)

Now I want to use the function "odeToVectorField".

F = odeToVectorField(a_dot_sym);

But I received the error.

Error using mupadengine/feval_internal

Unable to find the differential variables. Try specifying the dependent ODE variables using symbolic functions instead, such as 'syms

f(x)'.

Error in odeToVectorField>mupadOdeToVectorField (line 171)

T = feval_internal(symengine,'symobj::odeToVectorField',sys,x,stringInput);

Error in odeToVectorField (line 119)

sol = mupadOdeToVectorField(varargin);

Error in SED2 (line 82)

F = odeToVectorField(a_dot_sym);

What is the solution of this problem? Please help. Thanks in advance.

Answer 1

[Sam Chak]

I don't see any symbolic differential variables in these three so-called ODEs. That's why odeToVectorField() threw the error message. Can you rewrite the ODEs in terms of differential variables?

% -(100000*(4*a1 + 1/25)*(a1 - 99/100)^2 - 292000*a1*(a1 - 1)^2 + 292*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 - 800*a2*((53*a1)/500 - 499947/50000)^3*(a2 - 1)^2 - 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2 + 1200*a3*((53*a2*((53*a1)/5000 - 499947/500000))/500 + 10)^3*(a3 - 1)^2)*(292000*a1*(a1 - 1)^2 - 100000*(4*a1 + 1/25)*(a1 - 99/100)^2 - 292*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 + 800*a2*((53*a1)/500 - 499947/50000)^3*(a2 - 1)^2 + 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2 - 1200*a3*((53*a2*((53*a1)/5000 - 499947/500000))/500 + 10)^3*(a3 - 1)^2 + 1/100)
%  (108000*a1*(a1 - 1)^2 - 100*(8*a2 + 2/25)*((53*a1)/500 - 10)^3*(a2 - 99/100)^2 + 1200*a3*(a3 - 1)^2*((53*(2*a2 + 1/50)*((53*a1)/5000 - 1))/1000 + 10)^3 + 292*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 - 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2)*(108000*a1*(a1 - 1)^2 - 100*(8*a2 + 2/25)*((53*a1)/500 - 10)^3*(a2 - 99/100)^2 + 1200*a3*(a3 - 1)^2*((53*(2*a2 + 1/50)*((53*a1)/5000 - 1))/1000 + 10)^3 + 292*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 - 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2 - 1/100)
% -(108000*a1*(a1 - 1)^2 + 100*(12*a3 + 3/25)*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 99/100)^2 - 508*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 - 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2)*(508*a2*((53*a1)/500 - 10)^3*(a2 - 1)^2 - 100*(12*a3 + 3/25)*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 99/100)^2 - 108000*a1*(a1 - 1)^2 + 292*a3*((53*a2*((53*a1)/5000 - 1))/500 + 10)^3*(a3 - 1)^2 + 1/100)

A symbolic differential variable looks like this:

syms x(t)
x(t) = sin(t);
Dx   = diff(x,t)

Dx(t) = 

You can also find some examples in this documentation:

https://www.mathworks.com/help/symbolic/odetovectorfield.html

[Sol Elec]

This is my complete code.

% Define the constants
k = 0.73; 
c = [0.1, 0.1]; 
V1 = 10; 
N = 3; 

% Define the initial values and perturbation sizes 
a = [0.33, 0.33, 0.33];
da = [0.01, 0.01, 0.01]; 
Wd = [0.0053 0.0053 0.0053]; 

tau = 10; % Time constant
lambda = 0.01; % Safety margin

% Define the simulation time and step size
T = 100; % Simulation time
dt = 0.1; % Time step size

% Define a symbolic variable for time
syms t

% Define a symbolic vector 
a_sym = sym('a', [N, 1], 'real');

V_sym = sym('V', [N,1], 'real');

V_sym(1) = V1;

for j = 2:N
    V_sym(j) = V_sym(1)*(1-(1 - (1-2*a_sym(j-1) * (V_sym(j-1)/V_sym(1)) ))*Wd(j));
end

P_sym = sym('P', [N,1],'real');

for j = 1:N
    P_sym(j) = k*4*a_sym(j)*(1-a_sym(j))^2*V_sym(j)^3;
end


P_tot_sym = sum(P_sym);

dJ_da_sym = sym('dJ_da', [N, 1], 'real');

% Calculate the gradients of performance index using finite differences and symbolic expressions
for j = 1:N
    a_perturb_sym = a_sym;
    a_perturb_sym(j) = a_perturb_sym(j) + da(j);
    
    V_perturb_sym = sym('V_perturb', [N,1], 'real');
    V_perturb_sym(1) = V1;
    
    for k = 2:N
        V_perturb_sym(k) = V_perturb_sym(1)*(1-(1 - (1-2*a_perturb_sym(k-1) * (V_perturb_sym(k-1)/V_perturb_sym(1)) ))*Wd(k));
    end
    
    P_perturb_sym = sym('P_perturb', [N,1],'real');
    for k = 1:N
        P_perturb_sym(k) = k*4*a_perturb_sym(k)*(1-a_perturb_sym(k))^2*V_perturb_sym(k)^3;
    end
    
    P_tot_perturb_sym = sum(P_perturb_sym);
    
    dJ_da_sym(j) = (P_tot_perturb_sym - P_tot_sym)/da(j);
end

% rates of change 
a_dot_sym = sym('a_dot', [N, 1], 'real');

% Calculate the rates of change  using symbolic expressions
for j = 1:N
    a_dot_sym(j) = dJ_da_sym(j)*(dJ_da_sym(j) - lambda);
end

% Convert the second-order differential equation to a system of first-order differential equations using odeToVectorField function
F = odeToVectorField(a_dot_sym);

% Generate a MATLAB function from the system of first-order differential equations using matlabFunction function
f = matlabFunction(F,'vars', {'t','Y'});

% Define the initial conditions for the system of first-order differential equations
y0 = [a, zeros(1,N)];

% Solve the system of first-order differential equations using ode45 function
[t,y] = ode45(f,[0 T],y0);

% Plot the results
figure(1)
plot(t,y(:,1:N))
xlabel('Time (s)')
ylabel('value of a')
legend('a1','a2','a3')
title('Value of a vs time')

figure(2)
plot(t,sum(y(:,N+1:2*N),2))
xlabel('Time (s)')
ylabel('Total ')
title('Total vs time')

Please help

[Sam Chak]

Hi @Sol Elec

I couldn't find any symbolic differential variables in your entire code. Please try searching for the keyword 'diff('. I suggest that you directly display the three mathematical differential equations so that we can compare them with your code to ascertain the correctness of the ODEs in your code.