# Problem related with "odeToVectorField"

3 views (last 30 days)
Sol Elec on 2 Sep 2023
Commented: Sam Chak on 3 Sep 2023
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);
Unable to find the differential variables. Try specifying the dependent ODE variables using symbolic functions instead, such as 'syms
f(x)'.
T = feval_internal(symengine,'symobj::odeToVectorField',sys,x,stringInput);
Error in odeToVectorField (line 119)
Error in SED2 (line 82)
F = odeToVectorField(a_dot_sym);

Sam Chak on 2 Sep 2023
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:
Sol Elec on 3 Sep 2023
Moved: Sam Chak on 3 Sep 2023
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')
Sam Chak on 3 Sep 2023
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.