Runge Kutta Method for Matrix

24 visualizaciones (últimos 30 días)
Marcelo Boldt
Marcelo Boldt el 27 de En. de 2021
Comentada: Marcelo Boldt el 22 de Jun. de 2022
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
  1 comentario
Marcelo Boldt
Marcelo Boldt el 28 de En. de 2021
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

Iniciar sesión para comentar.

Respuestas (1)

Bharath Swaminathan
Bharath Swaminathan el 17 de Jun. de 2022
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.
  1 comentario
Marcelo Boldt
Marcelo Boldt el 22 de Jun. de 2022
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.

Iniciar sesión para comentar.

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Productos


Versión

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by