Runge Kutta Method for Matrix

16 次查看（过去 30 天）
Marcelo Boldt
Marcelo Boldt 2021-1-27
评论： Marcelo Boldt 2022-6-22
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 个评论
Marcelo Boldt
Marcelo Boldt 2021-1-28
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（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.
  1 个评论
Marcelo Boldt
Marcelo Boldt 2022-6-22
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Runge Kutta Methods 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by