This question is closed. Reopen it to edit or answer.

Vectorizing solve() for solving dynamical system for different parameter values

1 view (last 30 days)
atharva aalok
atharva aalok on 11 Sep 2021
Closed: atharva aalok on 8 Oct 2021
I will try to quickly explain the problem in brief and then pose the problem.
% Creating symbolic variables
% x0 y0 omega0 E0 are physical quantities
% Pm0 is a parameter generally in the range [0 1]
syms x0 y0 omega0 E0 Pm0
assume([x0 y0 omega0 E0], 'real');
assume(Pm0 ~= 0);
% x0 and y0 are related by
assumeAlso(x0*x0 + y0*y0 == 1);
% Equations for the dynamical system, (equating to 0 to find equilibrium
% points)
eqn1 = -y0*omega0 + x0*(1-x0*x0-y0*y0) == 0;
eqn2 = x0*omega0 + y0*(1-x0*x0-y0*y0) == 0;
eqn3 = -3.33*E0*y0 - 0.66*omega0 + 3.33 * Pm0 == 0;
eqn4 = 0.5*x0 - E0 + 0.5 == 0;
% solving the equations
eqn_final = [eqn1, eqn2, eqn3, eqn4];
variable_final = [x0, y0, omega0, E0];
a = solve(eqn_final, variable_final, "ReturnConditions", true);
% PmSol is a vector that contains different values for the parameter Pm0
PmSol = linspace(0.2, .5, 100);
% For each value of Pm0 from PmSol we need to solve the equations and find
% equilibrium pts.
% Then we need store the solutions for x0 to E0 for each Pm value.
% create vector of Fixed Point locations for each Pm val in PmSol
Fixed_points_locations = ones(2,length(PmSol)*4)*0;
% For each Pm value we will get either 2 solutions or no solutions
% if we get 2 solutions then we will store a vector like:
% [x0_1, y0_1 omega0_1 E0_1; x0_1, y0_1 omega0_1 E0_1] in above matrix. If
% no solutions then we will store zeros(2 by 4) in above matrix.
% THE AIM: the goal is to make the steps from here on out as fast as
% possible. This will become a part of a bigger problem needs to be
% lightning fast!
a = struct with fields:
x0: [1×1 sym] y0: [1×1 sym] omega0: [1×1 sym] E0: [1×1 sym] parameters: [1×1 sym] conditions: [1×1 sym]
ans = 
ans = 
Pm_val = PmSol;
% substitute all Pm values at once in a.conditions to get equation for
% parameter
subs_new = subs(a.conditions, Pm0, Pm_val);
ans = 
% Now all we need to do is solve each equation in subs_new (each of which
% looks like the above) get the parameter values, combine that with Pm_val
% and substitue in the expressions for x0 to E0 to get equilibrium points
% for all values of Pm
% THE ISSUE: HOW TO VECTORIZE THIS? (the following for loop)
for j = 1:length(PmSol)
% solve the each equation in subs_new for the parameter
parameter_val = solve(subs_new(j));
% make a vector of same dimensions as parameter vals with value
% Pm_val(j)
Pm_new = parameter_val * 0 + Pm_val(j);
% substitute parameter_val and Pm_new to get Equilibrium points
T1 = {[parameter_val], [Pm_new]};
T2 = {a.parameters, Pm0};
xFP = subs(a.x0, T2, T1);
yFP = subs(a.y0, T2, T1);
omegaFP = subs(a.omega0, T2, T1);
EFP = subs(a.E0, T2, T1);
[rows, cols] = size([xFP, omegaFP, EFP]);
% if no solutions add a zero matrix for the equilibrium point locations
if rows == 0
Fixed_points_locations(1:2,3*(j-1)+1:3*(j-1)+3) = ones(2, 3)*0;
Fixed_points_locations(1:2,3*(j-1)+1:3*(j-1)+3) = [xFP, omegaFP, EFP];

Answers (0)

Community Treasure Hunt

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

Start Hunting!