Runge Kutta Method for Matrix

Hello,

I have been developing a runge-kutta 4th order method to solve differential equations in matrix form (dx/dalpha = A x + Bu) where dapha stands for angle such that both A and B are functions of alpha. In addition, I have all the A matrices and B(alpha=0) = Identity matrix.

%% Laminate Conditions
ABD_Matrix_2 = -[227136.359975841,69646.0036239179,0,-454272.719951681,-139292.007247836,0;
    69646.0036239179,227136.359975841,7.27595761418343e-12,-139292.007247836,-454272.719951681,-1.63709046319127e-11;
    0,7.27595761418343e-12,78745.1781759614,-1.81898940354586e-12,-1.45519152283669e-11,-157490.356351923;
    -454272.719951681,-139292.007247836,0,1258428.76767918,370536.863822059,-7688.75830481177;
    -139292.007247836,-454272.719951681,-1.63709046319127e-11,370543.054681733,1166169.85888112,-7688.75830481172;
    -1.81898940354586e-12,-1.45519152283669e-11,-157490.356351923,-7687.72649486611,-7687.72649486608,419068.890196128];

%% Differentialmatrix for the dome of the tank (spherical geometry)
%% Creation of the Gesamtsystem
R = 2730;
s = pi/180; %step  

%% Insertion of ABD Matrix - Conditions

a = (ABD_Matrix_2(1,1)*ABD_Matrix_2(1,5) - ABD_Matrix_2(1,4)*ABD_Matrix_2(1,2));
a2 = (ABD_Matrix_2(1,1)*ABD_Matrix_2(1,5) - ABD_Matrix_2(1,4)*ABD_Matrix_2(1,2));
a4 = (- ABD_Matrix_2(2,1)*ABD_Matrix_2(4,4)*ABD_Matrix_2(1,2)+ ABD_Matrix_2(2,1)* ABD_Matrix_2(1,4)* ABD_Matrix_2(1,5)- ABD_Matrix_2(5,4)* ABD_Matrix_2(1,1)* ABD_Matrix_2(1,5)+ ABD_Matrix_2(2,4)*ABD_Matrix_2(1,4)*ABD_Matrix_2(1,2));
a11 = -ABD_Matrix_2(2,1)*ABD_Matrix_2(4,4) + ABD_Matrix_2(2,4)*ABD_Matrix_2(1,4);
a1 = - ABD_Matrix_2(1,2)*ABD_Matrix_2(4,4) + ABD_Matrix_2(2,4)*ABD_Matrix_2(1,4);
a12 = - ABD_Matrix_2(1,5)*ABD_Matrix_2(4,4) + ABD_Matrix_2(1,4)*ABD_Matrix_2(4,5);
a31 = - ABD_Matrix_2(1,2)* ABD_Matrix_2(1,4) + ABD_Matrix_2(2,4)* ABD_Matrix_2(1,1);
a33 = - ABD_Matrix_2(1,4)* ABD_Matrix_2(1,5) + ABD_Matrix_2(4,5)* ABD_Matrix_2(1,1);
det_dome = - ABD_Matrix_2(1,1)* ABD_Matrix_2(4,4) + ABD_Matrix_2(1,4)^2;
a41 = - ABD_Matrix_2(2,1)*ABD_Matrix_2(4,4)*ABD_Matrix_2(1,2) + ABD_Matrix_2(2,1)*ABD_Matrix_2(1,4)*ABD_Matrix_2(1,5) + ABD_Matrix_2(2,4)*ABD_Matrix_2(1,4)*ABD_Matrix_2(1,2) - ABD_Matrix_2(2,4)*ABD_Matrix_2(1,5)*ABD_Matrix_2(1,1);
a63 = - ABD_Matrix_2(2,4)*ABD_Matrix_2(4,4)*ABD_Matrix_2(1,5) + ABD_Matrix_2(4,5)*ABD_Matrix_2(2,4)*ABD_Matrix_2(1,4) + ABD_Matrix_2(5,4)*ABD_Matrix_2(1,4)*ABD_Matrix_2(1,5) - ABD_Matrix_2(5,4)*ABD_Matrix_2(1,1)*ABD_Matrix_2(4,5);
a46 = - ABD_Matrix_2(2,1)* ABD_Matrix_2(1,4) + ABD_Matrix_2(2,4)* ABD_Matrix_2(1,1);
a66 = - ABD_Matrix_2(2,4)* ABD_Matrix_2(1,4) + ABD_Matrix_2(5,4)* ABD_Matrix_2(1,1);
a64 = - ABD_Matrix_2(2,4)* ABD_Matrix_2(4,4) + ABD_Matrix_2(5,4)* ABD_Matrix_2(1,4);


for alpha = 1:89%L_sp
    phi = alpha*2/10;
    phi_rad = phi*pi/180; %angle in radians
    r =  2730 * cos(phi_rad);
    A_sp{alpha,:} =  [(a1 * sin(phi_rad))/(det_dome*r),1/2730 - (a1* cos(phi_rad))/(det_dome*r),(-sin(phi_rad)*a12)/(r*det_dome),-ABD_Matrix_2(4,4)/(det_dome*r),0,-ABD_Matrix_2(1,4)/(det_dome*r),0;
        -1/2730,0,1,0,0,0,0;
        (sin(phi_rad)*a31)/(r*det_dome),-(cos(phi_rad)*a31)/(r*det_dome),(-sin(phi_rad)*a33)/(r*det_dome),-ABD_Matrix_2(1,4)/(r*det_dome),0,-ABD_Matrix_2(1,1)/(r*det_dome),0;
        -sin(phi_rad)*(-ABD_Matrix_2(2,2)*sin(phi_rad)/r + sin(phi_rad)/det_dome*r*(a41)) ,- sin(phi_rad)*(ABD_Matrix_2(2,2)*cos(phi_rad)/r - cos(phi_rad)/det_dome*r*(a41)) , sin(phi_rad)*(ABD_Matrix_2(2,5)*-sin(phi_rad)/r + sin(phi_rad)/det_dome*r*(a41)) , ((-sin(phi_rad)*a11)/(det_dome*r)) , 1/2730 ,((sin(phi_rad)*a46)/(det_dome*r)),0;
        cos(phi_rad)*(ABD_Matrix_2(2,2)*-sin(phi_rad)/r + sin(phi_rad)/det_dome*r*(a41)) , cos(phi_rad)*(ABD_Matrix_2(2,2)*cos(phi_rad)/r - cos(phi_rad)/det_dome*r*(a41)), cos(phi_rad)*(ABD_Matrix_2(2,5)*sin(phi_rad)/r - sin(phi_rad)/det_dome*r*(a41)), -1/2730 + (cos(phi_rad) * a11)/(det_dome*r) ,0,-((cos(phi_rad)*a46)/(det_dome*r)),-0.4*(0.5-0.5*0.009);
        sin(phi_rad)*(ABD_Matrix_2(2,5)*-sin(phi_rad)/r + sin(phi_rad)/det_dome*r*(a41)) , sin(phi_rad)*(ABD_Matrix_2(2,5)*cos(phi_rad)/r - cos(phi_rad)/det_dome*r*(a41)), sin(phi_rad)*(sin(phi_rad)*(ABD_Matrix_2(5,5)*r)-(sin(phi_rad)/(det_dome*r))*a63) + 546 * r, ((sin(phi_rad)*a64)/(det_dome*r)), -1, ((-sin(phi_rad)*a66)/(det_dome*r)) ,0;
        0,0,0,0,0,0,0];
    
    
end


Differentialmatrix_Kugel_final = cell2mat(A_sp);


%% Runge Kutta 4th Order - Tank - Spherical Dome
B{1,:} = eye(7,7)
for j= 1:88-1 %             % calculation loop
    %j = j +25;
    m_1 = A_sp{j,:}*B_sp{j,:}  ;                             % calculating coefficient
    m_2 = (A_sp{j,:}+A_sp{j+1,:}/2)*(B_sp{j,:}+0.5*0.2*m_1); % for replacement
    m_3 = (A_sp{j,:}+A_sp{j+1,:}/2)*(B_sp{j,:}+0.5*0.2*m_2);
    m_4 = (A_sp{j,:})*(B_sp{j,:}+m_3*0.2);
    B_sp{j+1,:} = B_sp{j,:} + (0.2/6)*(m_1+2*m_2+2*m_3+m_4);  % main equation
%     Ubertragungsmatrixcell_sp{j,:} = B_sp{:,j};
    
end

Unfortunately my integration does not converge and I dont know why, do you have any suggestion regarding this problem?

Thank you

[Marcelo Boldt]

In an attempt of improving the algorithm, I re wrote it as:

B_sp = cell(89,1);
B_sp{1,:} = eye(7,7);
TMMdot = [B_sp,A_sp]; 
B{1,:} = B_sp{1,:};

for j= 1:88-1 %      712-1                     %25:1:89-1        % calculation loop
    %j = j +25;
    m_1 = TMMdot(B{j,:},A_sp{j,:});                             % calculating coefficient
    m_2 = TMMdot( B{j,:} + 0.1 *  m_1,A_sp{j+1,:});
    m_3 = TMMdot( B{j,:} + 0.1 *  m_1,A_sp{j+1,:});
    m_4 = TMMdot( B{j,:} + 0.1 *  m_1,A_sp{j+2,:});
    B{j+1,:} = B{j,:} + (0.2/6)*(m_1+2*m_2+2*m_3+m_4);  % main equation
%     Ubertragungsmatrixcell_sp{j,:} = B_sp{:,j};
    
end

But unfortunately I am getting the following mistake:

Index in position 1 is invalid. Array indices must be positive integers or logical values.

Help will be highly appreciated

Answer 1

[Bharath Swaminathan]

TMMdot is an array consisting of B_sp and A_sp. You can only pass integers as indices, but you are passing B{j,:} and A_sp{j,:} which is incorrect.

Having said that, i haven't read your problem completely. If you want to solve differential equations, you can use Matlab's ode45, ode15s, ode23s etc. solvers. If you are trying to implement a Runge-Kutta solver from scratch, then there are lot of online resources which give the correct implementation for the RK4 method. Your implementations has a lot of flaws - why are you multiplying A_sp{j,:} and B_sp{j,:}? the equation is xdot = A_sp.x +B_sp.u right? Your expressions for m1,m2,m3,m4 needs to be revised.

[Marcelo Boldt]

Hi, Thank you for your answer. I will go through it and see how it goes. I agree with you that some expressions are not solid, or either a false, however these expressions have been given with a high level of accuracy. So wisest step is to revise them and check for possible mistakes.