Contents
function out = LQR_LTV_smallscale_splitting(method, istest)
if nargin < 1
method.order = 2;
method.additive = false;
method.symmetric = false;
end
if nargin < 2
istest = 0;
end
Nx = 10;
Nx2 = Nx^2;
xx = linspace(0,1, Nx+2);
x = xx(2:end-1);
dx = 1/(Nx+1);
e = ones(Nx, 1);
A_1D = 1/dx^2*spdiags([e, -2*e, e], -1:1, Nx, Nx);
[X, Y] = meshgrid(x,x);
I_1D = eye(Nx, Nx);
A_2D = kron(A_1D, I_1D) + kron(I_1D, A_1D);
M = speye(Nx2, Nx2);
Nb = 3;
B = zeros(Nx2, Nb);
for j = 1:Nb
B(:, j) = reshape((X > j/4) & (X < j/4 + 1/8) & (Y > j/4) & (Y < j/4 + 1/8), Nx2, 1);
end
C = 1/Nx2 * ones(1, Nx2);
alpha = @(t) 1 + 10*sin(2*pi*t);
mu = @(t) 2 + 7.1*sin(2*pi*t);
dmu = @(t) 7.1*2*pi*cos(2*pi*t);
beta = @(t) 3 + cos(t);
gamma = @(t) 1 - min(t, 1);
At = @(t) alpha(t) * A_2D;
Mt = @(t) mu(t) * M;
dMt = @(t) dmu(t) * M;
Bt = @(t) beta(t) * B;
Ct = @(t) gamma(t) * C;
eqn.A_time = At;
eqn.E_time = Mt;
eqn.dt_E_time = dMt;
eqn.B_time = Bt;
eqn.C_time = Ct;
eqn.haveE = 1;
eqn.type = 'T';
L0 = zeros(Nx2, 1);
eqn.L0 = L0;
eqn.D0 = eye(size(L0, 2));
eqn.Rinv = 1;
oper = operatormanager('default');
eqn.LTV = 1;
t0 = 0;
tend = 0.1;
Nt = 100;
General splitting parameters
opts.splitting.time_steps = linspace(t0, tend, Nt+1);
opts.splitting.order = method.order;
opts.splitting.additive = method.additive;
opts.splitting.symmetric = method.symmetric;
opts.splitting.info = 2;
opts.splitting.intermediates = 1;
opts.splitting.trunc_tol = eps;
opts.splitting.quadrature.type = 'adaptive';
opts.splitting.quadrature.tol=1e-8 ;
Matrix exponential actions
opts.exp_action.method = 'LTV';
opts.exp_action.tol = 1e-8;
Compute the approximation
tic;
[out, ~,opts, ~] = mess_splitting_dre(eqn,opts,oper);
toc;
Step: 1 of 100, time = 0
Step: 2 of 100, time = 1.000000e-03
Step: 3 of 100, time = 2.000000e-03
Step: 4 of 100, time = 3.000000e-03
Step: 5 of 100, time = 4.000000e-03
Step: 6 of 100, time = 5.000000e-03
Step: 7 of 100, time = 6.000000e-03
Step: 8 of 100, time = 7.000000e-03
Step: 9 of 100, time = 8.000000e-03
Step: 10 of 100, time = 9.000000e-03
Step: 11 of 100, time = 1.000000e-02
Step: 12 of 100, time = 1.100000e-02
Step: 13 of 100, time = 1.200000e-02
Step: 14 of 100, time = 1.300000e-02
Step: 15 of 100, time = 1.400000e-02
Step: 16 of 100, time = 1.500000e-02
Step: 17 of 100, time = 1.600000e-02
Step: 18 of 100, time = 1.700000e-02
Step: 19 of 100, time = 1.800000e-02
Step: 20 of 100, time = 1.900000e-02
Step: 21 of 100, time = 2.000000e-02
Step: 22 of 100, time = 2.100000e-02
Step: 23 of 100, time = 2.200000e-02
Step: 24 of 100, time = 2.300000e-02
Step: 25 of 100, time = 2.400000e-02
Step: 26 of 100, time = 2.500000e-02
Step: 27 of 100, time = 2.600000e-02
Step: 28 of 100, time = 2.700000e-02
Step: 29 of 100, time = 2.800000e-02
Step: 30 of 100, time = 2.900000e-02
Step: 31 of 100, time = 3.000000e-02
Step: 32 of 100, time = 3.100000e-02
Step: 33 of 100, time = 3.200000e-02
Step: 34 of 100, time = 3.300000e-02
Step: 35 of 100, time = 3.400000e-02
Step: 36 of 100, time = 3.500000e-02
Step: 37 of 100, time = 3.600000e-02
Step: 38 of 100, time = 3.700000e-02
Step: 39 of 100, time = 3.800000e-02
Step: 40 of 100, time = 3.900000e-02
Step: 41 of 100, time = 4.000000e-02
Step: 42 of 100, time = 4.100000e-02
Step: 43 of 100, time = 4.200000e-02
Step: 44 of 100, time = 4.300000e-02
Step: 45 of 100, time = 4.400000e-02
Step: 46 of 100, time = 4.500000e-02
Step: 47 of 100, time = 4.600000e-02
Step: 48 of 100, time = 4.700000e-02
Step: 49 of 100, time = 4.800000e-02
Step: 50 of 100, time = 4.900000e-02
Step: 51 of 100, time = 5.000000e-02
Step: 52 of 100, time = 5.100000e-02
Step: 53 of 100, time = 5.200000e-02
Step: 54 of 100, time = 5.300000e-02
Step: 55 of 100, time = 5.400000e-02
Step: 56 of 100, time = 5.500000e-02
Step: 57 of 100, time = 5.600000e-02
Step: 58 of 100, time = 5.700000e-02
Step: 59 of 100, time = 5.800000e-02
Step: 60 of 100, time = 5.900000e-02
Step: 61 of 100, time = 6.000000e-02
Step: 62 of 100, time = 6.100000e-02
Step: 63 of 100, time = 6.200000e-02
Step: 64 of 100, time = 6.300000e-02
Step: 65 of 100, time = 6.400000e-02
Step: 66 of 100, time = 6.500000e-02
Step: 67 of 100, time = 6.600000e-02
Step: 68 of 100, time = 6.700000e-02
Step: 69 of 100, time = 6.800000e-02
Step: 70 of 100, time = 6.900000e-02
Step: 71 of 100, time = 7.000000e-02
Step: 72 of 100, time = 7.100000e-02
Step: 73 of 100, time = 7.200000e-02
Step: 74 of 100, time = 7.300000e-02
Step: 75 of 100, time = 7.400000e-02
Step: 76 of 100, time = 7.500000e-02
Step: 77 of 100, time = 7.600000e-02
Step: 78 of 100, time = 7.700000e-02
Step: 79 of 100, time = 7.800000e-02
Step: 80 of 100, time = 7.900000e-02
Step: 81 of 100, time = 8.000000e-02
Step: 82 of 100, time = 8.100000e-02
Step: 83 of 100, time = 8.200000e-02
Step: 84 of 100, time = 8.300000e-02
Step: 85 of 100, time = 8.400000e-02
Step: 86 of 100, time = 8.500000e-02
Step: 87 of 100, time = 8.600000e-02
Step: 88 of 100, time = 8.700000e-02
Step: 89 of 100, time = 8.800000e-02
Step: 90 of 100, time = 8.900000e-02
Step: 91 of 100, time = 9.000000e-02
Step: 92 of 100, time = 9.100000e-02
Step: 93 of 100, time = 9.200000e-02
Step: 94 of 100, time = 9.300000e-02
Step: 95 of 100, time = 9.400000e-02
Step: 96 of 100, time = 9.500000e-02
Step: 97 of 100, time = 9.600000e-02
Step: 98 of 100, time = 9.700000e-02
Step: 99 of 100, time = 9.800000e-02
Step: 100 of 100, time = 9.900000e-02
Elapsed time is 6.927538 seconds.
if not(istest)
t = opts.splitting.time_steps;
figure;
plot(t, out.ms);
title('Ranks of approximations over time');
y = zeros(1,length(out.Ks));
for i=1:length(out.Ks)
y(i) = out.Ks{i}(1,1);
end
figure;
plot(t, y);
title('evolution of component (1,1) of the optimal feedback');
else
if abs(norm(out.Ds{1}) / 6.078945091766749e-05 - 1) >= 1e-10
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
end
ans =
struct with fields:
Ls: {1×101 cell}
Ds: {1×101 cell}
ms: [101×1 double]
Ks: {1×101 cell}