Problem related with "odeToVectorField"

7 次查看（过去 30 天）
Sol Elec
Sol Elec 2023-9-2
评论： Sam Chak 2023-9-3
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.

请先登录，再回答此问题。

回答（1 个）

Sam Chak
Sam Chak 2023-9-2
Ran in:
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
  2 个评论
Sol Elec
Sol Elec 2023-9-3
移动：Sam Chak 2023-9-3
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
Sam Chak 2023-9-3
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Equation Solving 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by