Метод интеграции

Question

Я занимаюсь исследованием, которое требует точности вычисления и размера шага, в котором он необходим. У меня есть собственные значения систем, но я не знаю, как рассчитать размер шага по собственным значениям точности. Я довел шаг по времени до минимума, чтобы он мог представлять собой аналитическое решение (y2), но, поскольку размеры шагов различны, я не могу определить точность. Любой имеет какие-либо идеи или может помочь мне найти код Matlab, реализующий эти вещи.

clear all

Rs = 1;

E = 5;
Rc = 10;
C1 = 0.09; C2 = 0.09; C3 = 0.09; C4 = 0.09; C5 = 0.09; C6 = 0.09; C7 = 0.09; C8 = 0.01; C9 = 0.01; C10 = 0.9;
L1 = 0.001; L2 = 0.001; L3= 0.001; L4 = 0.0001; L5 = 0.001; L6 = 0.001; L7 = 0.001; L8 = 0.001; L9 = 0.001; L10 = 0.01;

h = 0.0003;
A=[-Rs/L1,0,0,0,0,0,0,0,0,0,-1/L1,0,0,0,0,0,0,0,0,0;
    0,-Rs/L2,0,0,0,0,0,0,0,0, 1/L2,-1/L2,0,0,0,0,0,0,0,0;
    0,0,-Rs/L3,0,0,0,0,0,0,0,0,1/L3,-1/L3,0,0,0,0,0,0,0;
    0,0,0,-Rs/L4,0,0,0,0,0,0,0,0,1/L4,-1/L4,0,0,0,0,0,0;
    0,0,0,0,-Rs/L5,0,0,0,0,0,0,0,0,1/L5,-1/L5,0,0,0,0,0;
    0,0,0,0,0,-Rs/L6,0,0,0,0,0,0,0,0,1/L6,-1/L6,0,0,0,0;
    0,0,0,0,0,0,-Rs/L7,0,0,0,0,0,0,0,0,1/L7,-1/L7,0,0,0;
    0,0,0,0,0,0,0,-Rs/L8,0,0,0,0,0,0,0,0,1/L8,-1/L8,0,0;
    0,0,0,0,0,0,0,0,-Rs/L9,0,0,0,0,0,0,0,0,1/L9,-1/L9,0;
    0,0,0,0,0,0,0,0,0,-Rs/L10,0,0,0,0,0,0,0,0,1/L10,-1/L10;
    1/C1,-1/C1,0,0,0,0,0,0,0,0,-1/(Rc*C1),0,0,0,0,0,0,0,0,0;
    0,1/C2,-1/C2,0,0,0,0,0,0,0,0,-1/(Rc*C2),0,0,0,0,0,0,0,0;
    0,0,1/C3,-1/C3,0,0,0,0,0,0,0,0,-1/(Rc*C3),0,0,0,0,0,0,0;
    0,0,0,1/C4,-1/C4,0,0,0,0,0,0,0,0,-1/(Rc*C4),0,0,0,0,0,0;
    0,0,0,0,1/C5,-1/C5,0,0,0,0,0,0,0,0,-1/(Rc*C5),0,0,0,0,0;
    0,0,0,0,0,1/C6,-1/C6,0,0,0,0,0,0,0,0,-1/(Rc*C6),0,0,0,0;
    0,0,0,0,0,0,1/C7,-1/C7,0,0,0,0,0,0,0,0,-1/(Rc*C7),0,0,0;
    0,0,0,0,0,0,0,1/C8,-1/C8,0,0,0,0,0,0,0,0,-1/(Rc*C8),0,0;
    0,0,0,0,0,0,0,0,1/C9,-1/C9,0,0,0,0,0,0,0,0,-1/(Rc*C9),0;
    0,0,0,0,0,0,0,0,0,1/C10,0,0,0,0,0,0,0,0,0,-1/(Rc*C10)];

B = [   1/L1 ,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0   ]'; 
u = 5;

Apr = inv(eye(20) -A*h);
Bpr = Apr* h * B ;

tstop = 40.0;

time = 0:h:tstop;
x = zeros(20, length(time));
tic
for k=2:length(time)
    if time(k)>8
        u=10;
    end
    x(:,k) = Apr * x(:,k-1) + Bpr * u;

end
toc

h2=0.000002;
Apr2 = inv(eye(20) -A*h2);
Bpr2 = Apr2* h2 * B ;
time2 = 0:h2:tstop;
y2 = zeros(20, length(time2));
tic
u2=5;
for z=2:length(time2)
    if time2(z)>8
        u2=10;
    end
    y2(:,z) = Apr2 * y2(:,z-1) + Bpr2 * u2;

end
toc
% -----------------------------------------
figure(1)
e=eig(A);
r=max(e)/min(e);
circle=e*h;
plot(real(circle), imag(circle), '*m' ,'MarkerSize', 10); hold on;
plot(real(e), imag(e), '*r' ,'MarkerSize', 10); hold on;
figure(2)
plot(time, x(20,:), 'm--');hold on;
plot(time2, y2(20,:), 'g--');