Tartalomjegyzék

Script d2017_09_11_maklc_v2_unc

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.

Output:

┌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

preprocessing, analysis

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);

Output:

│   ┌P_generate_symvars_v5
│   └ 0.4835 [sec]

THE ALGORITHM

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;

Output:

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)

Output:

│   ┌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

Simulation

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.

Published with MATLAB® R2015a, and sanitarized by BeautifulSoup (Python script using BSoup is written by Polcz).