TELJES MATLAB SCRIPT KIEGÉSZÍTŐ FÜGGVÉNYEKKEL
file: d2017_09_11_maklc_v2_unc.m author: Peter Polcz <ppolcz@gmail.com>
Created on 2017. September 11.
┌d2017_09_11_maklc_v3_unc │ - Persistence for `d2017_09_11_maklc_v3_unc` reused (inherited) [run ID: 8332, 2017.09.18. Monday, 09:42:54] │ - Script `d2017_09_11_maklc_v3_unc` backuped
P_generate_symvars_v5(2,1)
syms v1 k2 k3 real
xeq = [
k3^2/(k2*v1)
v1/k3
];
fbar_sym = simplify([
v1 - k2*x1*x2^2
k2*x1*x2^2 - k3*x2
]);
f_sym = subs(fbar_sym, x, x+xeq);
f_sym = simplify(expand(f_sym));
f_sym = collect(f_sym, x);
│ ┌P_generate_symvars_v5 │ └ 0.4835 [sec]
n = 2;
v1 = 4;
k2 = d1;
k3 = 2;
[b_sym,a_sym] = numden(f_sym);
fxo = f_sym;
A_sym = sym(zeros(n,n));
for ii = 1:n
for jj = 1:n
[coefficients,bases] = coeffs(b_sym(ii), x(jj));
% which elements of coefficients contains some power of x(jj)
% - subs(t, x(jj), 0) where x(jj) exists they will be 0
% - sign(abs(...)) where nonzero will be 1
tmp = 1 - sign(abs(subs(bases, x(jj), 0)));
% add those elements to A(ii,jj) where check is true, divided by x(jj)
A_sym(ii,jj) = sum(tmp .* coefficients .* (bases/x(jj))) / a_sym(ii);
% remove these elements from f and b
f_sym(ii) = simplify(f_sym(ii) - A_sym(ii,jj)*x(jj));
b_sym(ii) = simplify(b_sym(ii) - A_sym(ii,jj)*a_sym(ii)*x(jj));
end
end
f_sym = fxo;
clear fxo tmp coefficients bases
pvpa('A(x)x - f(x)', simplify(A_sym * x - f_sym))
pcz_display(A_sym)
fbar_sym = subs(fbar_sym);
f_sym = subs(f_sym);
A_sym = subs(A_sym);
xeq = subs(xeq);
pcz_display(A_sym)
sdps = sdpsettings('solver','mosek', 'debug',false, 'verbose',true);
sdps.mosek.MSK_IPAR_NUM_THREADS = 8;
sdps.alphamethod = false;
A(x)x - f(x) [size=2x1,disp] =
0
0
A_sym [2x2] =
[ -(k2*k3^2*v1*x2^2 + 2*k2*k3*v1^2*x2 + k2*v1^3)/(k3^2*v1), -(x2*k3^4 + 2*v1*k3^3)/(k3^2*v1)]
[ (k2*k3^2*v1*x2^2 + 2*k2*k3*v1^2*x2 + k2*v1^3)/(k3^2*v1), (x2*k3^4 + v1*k3^3)/(k3^2*v1)]
A_sym [2x2] =
[ - d1*x2^2 - 4*d1*x2 - 4*d1, - x2 - 4]
[ d1*x2^2 + 4*d1*x2 + 4*d1, x2 + 2]
% -1: just sym2lfr
% 0: minlfr
% 1: v5 - split/merge
% 2: v6 - LFR reduction
method = 0;
A_sym = simplify(A_sym);
A_lfr = sym2lfr(A_sym);
model = v2struct(A_lfr, f_sym, x_n, d_n, Dd_n);
model.vars = {x,d,Dd};
if method == 0
model = P_algo_part1_v5_testing(model, @P_generate_matrices_v5);
elseif method == 1
model = P_algo_part1_v5(model, @P_generate_matrices_v5);
elseif method == 2
model = P_algo_part1_v6(model, @P_generate_matrices_v5);
end
expand(model.Pi)
% pcz_num2str_latex(model.lfr.M22)
% pcz_num2str_latex(model.lfr.M21)
% pcz_num2str_latex(model.lfr.M12)
% pcz_num2str_latex(model.lfr.M11)
% pcz_num2str_latex(model.lfr.A)
% pcz_num2str_latex(model.lfr.B)
% pcz_num2str_latex(model.lfr.C)
% pcz_num2str_latex(model.lfr.D)
│ ┌P_algo_part1_v5_testing
│ │ - Persistence for `P_algo_part1_v5_testing` initialized [run ID: 3113, 2017.09.18. Monday, 09:51:49]
│ │ - Script `P_algo_part1_v5_testing` backuped
│ │ ┌P_model2T16_v5
│ │ │ ┌P_lfrdata_v1
│ │ │ └ 0.010632 [sec]
│ │ └ 0.31216 [sec]
│ │ ┌P_generate_matrices_v5
│ │ │ - Persistence for `P_generate_matrices_v5` initialized [run ID: 4178, 2017.09.18. Monday, 09:51:50]
│ │ │ - Script `P_generate_matrices_v5` backuped
│ │ │ ┌P_annihilator_linsolve_v5
│ │ │ │ - Persistence for `P_annihilator_linsolve_v5` initialized [run ID: 8390, 2017.09.18. Monday, 09:51:50]
│ │ │ │ - Script `P_annihilator_linsolve_v5` backuped
│ │ │ │ ┌P_annihilator_linsolve_v5 - core computation
│ │ │ │ └ 0.15626 [sec]
│ │ │ │ ┌P_annihilator_linsolve_v5 - quick check I. - obtained annihilator, prec: 10
│ │ │ │ │ - maximal difference: 1.580225e-17, row nr. 1
│ │ │ │ └ 0.10644 [sec]
│ │ │ │ ┌P_annihilator_linsolve_v5 - beautify annihilator + !roundings!
│ │ │ │ └ 0.009491 [sec]
│ │ │ │ ┌P_annihilator_linsolve_v5 - quick check II. - the beautified annihilator, prec: 10
│ │ │ │ │ - maximal difference: 8.104628e-15, row nr. 3
│ │ │ │ └ 0.097852 [sec]
│ │ │ └ 0.38448 [sec]
│ │ │ ┌P_row_reduction
│ │ │ └ 0.33576 [sec]
│ │ │ [FAILED] M is a true annihilator of Pia `simplify(M * Pia)`
│ │ │ [ OK ] M is practically a true annihilator of Pia [using a random numerical substitution (<1e-10)]
│ │ └ 3.0333 [sec]
│ └ 3.3668 [sec]
ans =
-2^(1/2)*d1*x1
(1122^(1/2)*d1*x1*x2)/1122 - (4*1122^(1/2)*x2^2)/561 - (4*1122^(1/2)*d1*x1*x2^2)/561
(66^(1/2)*x2^2)/33 + (4*66^(1/2)*d1*x1*x2)/33 + (66^(1/2)*d1*x1*x2^2)/33
return
X_chose = 'X0';
X.X0 = [
-0.833 1.4
-0.7 1.3
];
X.X1 = [
-1.4 1.4
-0.7 1.3
];
X.X2 = [
-1.9 1.9
-0.7 1.5
];
d_lims = [
0.8 1.2
];
res = 300;
resd = 30;
d_span = linspace(d_lims(1), d_lims(2), resd);
X = X.(X_chose);
model.lims = { X d_lims 0 };
xeq_fh = matlabFunction(xeq + eps*d1);
xeq_num = xeq_fh(d_span);
xeq_min = min(xeq_num,[],2);
xeq_max = max(xeq_num,[],2);
margin = 0.2;
X_minmax = [ X(:,1) + xeq_min - margin, X(:,2) + xeq_max + margin ];
[domains, result] = P_ndsolve_v5(model,sdps);
P = double(result.P);
V = model.Pib' * P * model.Pib;
V = subs(V,x,x-xeq);
V_fh = matlabFunction(V, 'vars', {x1 x2 d1});
[x1_num,x2_num,d_num] = ndgrid(...
linspace(0,X1(1,2)+1.2,res),...
linspace(X1(2,1)-0.2,X1(1,2)+0.2,res) + 2,...
d_span);
V_num = V_fh(x1_num, x2_num, d_num);
Vi = max(V_num,[],3);
Vo = min(V_num,[],3);
C = pcontour(x1_num(:,:,1), x2_num(:,:,1), Vi, [1,1]); grid on; hold on;
Co = pcontour(x1_num(:,:,1), x2_num(:,:,1), Vo, [1,1]); grid on; hold on;
I = [3 4 2 1 3];
% plot(domains.X_v(I,1), domains.X_v(I,2), '--ob')
disp 'inner'
for i = 1:numel(C)
C(i).area = polyarea(C(i).contour(1,:), C(i).contour(2,:));
disp(C(i).area)
end
disp 'outer'
for i = 1:numel(Co)
Co(i).area = polyarea(Co(i).contour(1,:), Co(i).contour(2,:));
disp(Co(i).area)
end
xeq_eps = xeq;
xeq_eps(2) = xeq(2) + eps*d1;
xeq_fh = matlabFunction(xeq_eps);
xeq_num = xeq_fh(linspace(d_lims(1), d_lims(2), resd));
plot(xeq_num(1,:),xeq_num(2,:),'b.')
result.sol
result.nr_Variables
numel(model.Pib)
numel(model.Pia)
return
fbar_fh = matlabFunction(subs(fbar_sym,d,d_num), 'vars', {'t', x});
f_fh = matlabFunction(subs(f_sym,d,d_num), 'vars', {'t', x});
% figure
% [tt,xx] = ode45(f_fh, [0 10], [-1.1,-0.7]);
% plot(xx(:,1),xx(:,2),'-')
figure, hold on
for d_num = d_lims(1):0.1:d_lims(2)
fbar_fh = matlabFunction(subs(fbar_sym,d,d_num), 'vars', {'t', x});
[tt,xx] = ode45(fbar_fh, [0 10], [-0.1,1.5]);
plot(xx(:,1),xx(:,2),'.-')
plot( [ 0 0 max(xx(:,1)) ] , [ max(xx(:,2)) 0 0 ] , 'k', 'LineWidth', 3)
end
End of the script.