Trouble formulating Range Kutta

Hello there, I'm having trouble coming with a way to solve this resolving range Kutta, to give me a graphic using the commented part in the end of the code. I am trying to solve this Matrix equation system of third order.

with this code

clear all, close all, clc
s = 7.5;                    % semi span
c = 2;                      % chord
hline = 80;                 % hinge line as percentage of chord
xh = hline / 100 * c;       % distance of hinge line aft of leading edge
fa = 40;                    % elastic axis as percentage of chord
xf = fa / 100 * c;          % position of elastic axis aft of le
xac = 0.25 * c;             % distance of ac aft of le
xe = xf - xac;              % distance of elastic axis aft of ac
e = xe / c;                 % eccentricity of fa aft of ac per chord
heave = 3;                  % heave angle
alpha = 3;                  % pitch angle
beta = 3;                   % control surface angle


chat = 2 * (xh / c) - 1;

T4 = - acos(chat) + chat * sqrt(1 - chat ^ 2);
T10 = sqrt(1 - chat ^ 2) + acos(chat);
T11 = acos(chat) * (1 - 2 * chat) + sqrt(1 - chat ^ 2) * (2 - chat);
T12 = sqrt(1 - chat ^ 2) * (2 + chat) - acos(chat) * (2 * chat + 1);

% Aerodynamic coefficients

aw = 2 * pi;                % lift per incidence
ac = aw * T10 / pi;         % lift per control rotation
bw = e * aw;                % pitching moment per control rotation
bc = e * aw * T10 / pi;     % pitching moment per control rotation
cw = - T12 / 2;             % hinge moment per incidence
cc = - T12 * T10 / 2 / pi;  % hinge moment per control rotation

Mtdot = -1.2;               % unsteady torsional damping term
Mbdot = -0.1;               % unsteady control rotation damping term

rho = 1.225;                % air density

% Mass data
m = 400;                    % mass per unit area
mmain = m;                  % mass per unit area of main surface
mcont = m;                  % mass per unit area of control surface

% Stiffness data
EI =  4e+7;                 % flexural rigidity
GJ =  8e+6;                 % torsional rigidity
kbeta = 1e+4;               % control stiffness per unit span
kh = 4 * EI / s^3;          % bending spring stiffness
ka = GJ / s;                % torsion spring stiffness
ka3 = 100*ka*alpha^3;       %non linear stiffness

GJstr = sprintf('%0.5g',GJ); EIstr = sprintf('%0.5g',EI);
astring = ['xf/c = ', num2str(xf/c),' xh/c = ',num2str(xh/c),' EI = ',EIstr,...
    ' GJ = ',GJstr,' kbeta = ', num2str(kbeta)];
disp(astring)

modes = 3;

% Inertia matrix
A = zeros(modes, modes);
% Main surface
Amain = zeros(modes, modes);
Amain(1,1) = xh * s / 5;
Amain(1,2) = (xh^2 / 2 - xh * xf) * s / 4;
Amain(2,1) = Amain(1,2);
Amain(2,2) = (xh^3 / 3 - xh^2 * xf + xf^2 * xh) * s / 3;
% Control surface
Acont = zeros(modes, modes);
Acont(1,1) = (c - xh) * s / 5;
Acont(1,2) = ((c^2 - xh^2) / 2 - xf * (c - xh)) * s / 4;
Acont(1,3) = ((c^2 - xh^2) / 2 - xh * (c - xh)) * s / 3;
Acont(2,2) = ((c^3 - xh^3) / 3 - xf * (c^2 - xh^2) + xf^2 * (c - xh)) * s / 3;
Acont(2,3) = ((c^3 - xh^3) / 3 - xf * (c^2 - xh^2) / 2 - xh * (c^2 - xh^2) / 2 + xf * xh * (c - xh)) * s / 2;
Acont(3,3) = ((c^3 - xh^3) / 3 - xh * (c^2 - xh^2) + xh^2 * (c - xh)) * s;
Acont(2,1) = Acont(1,2);
Acont(3,1) = Acont(1,3);
Acont(3,2) = Acont(2,3);

A = mmain * Amain + mcont * Acont;

% Structural damping matrix is zero
D = zeros(modes, modes);

B = zeros(modes, modes);     % aero damping - based on rho*V*B
C = zeros(modes, modes);     % aero stiffness - based on rho*V^2*C

B(1,1) = aw * c * s / 10;
B(1,2) = 0;
B(1,3) = 0;
B(2,1) = - bw * c^2 * s / 8;
B(2,2) = - Mtdot * c^3 * s / 24;
B(2,3) = 0;
B(3,1) = - cw * c^2 * s / 6;
B(3,2) = 0;
B(3,3) = - Mbdot * c^3 * s / 8;

C(1,1) = 0;
C(1,2) = aw * c * s / 8;
C(1,3) = ac * c * s / 6;
C(2,1) = 0;
C(2,2) = - bw * c^2 * s / 6;
C(2,3) = - bc * c^2 * s / 4;
C(3,1) = 0;
C(3,2) = - cw * c^2 * s / 4;
C(3,3) = - cc * c^2 * s / 2;


E = zeros(modes, modes);
E(1,1) = kh; 
E(2,2) = ka;
E(3,3) = kbeta * s;
n = 3;
V = 117;
t = linspace(0,10,500);
Mat=[zeros(n,n),eye(n,n);
    -A\(rho*V*V*C+E), -A\(rho*V*B+D)];

X1 = [x(1)
    x(2)
    x(3)
    x(4)
    x(5)
    x(6)];

x(1) = x1;
x(2) = x2;
x(3) = x3;
x(4) = x4;
x(5) = x5;
x(6) = x6;

%não linearidade em alpha

X2 = [0
    x(2)^3
    0
    0
    0
    0];



Mat3 = @(t,x) Mat.*X1 + [-ka3*A\[0; 1 ;0] zeros(3,1)]*X2;

% a=Mat;
% 
% b=[];
% 
% c=eye(length(Mat));
% 
% d=[];
% 
% sys=ss(a,b,c,d);
% 
% x0=[10e-2 0 0 0 0 0];
% 
% [y,t,x]=initial(sys,x0,t);
% 
% figure(1)
% 
% subplot(2,1,1)
% plot(t,y(:,1),'b',t,y(:,2),'r',t,y(:,3),'y')
% grid on
% xlabel('Tempo(s)')
% ylabel('Deslocamento[m]')
% % title(str)
% % h1=legend('$\kappa- 90m/s$ ','$\theta- 90m/s$','$\kappa- 30m/s$ ','$\theta- 30m/s$');
% % set(h1, 'Interpreter', 'latex');
% subplot(2,1,2)
% plot(t,y(:,3),'b',t,y(:,4),'r',t,y(:,6),'y')
% grid on
% xlabel('Tempo(s)')
% ylabel('Velocidade[m/s]')
% hold on;

Thanks in advance!

Answer 1

[Torsten]

The expression

[-ka3*A\[0; 1 ;0] zeros(3,1)]

must be 6x6 instead of 3x2 in order to work for your system.

And the numerical integration method is "Runge-Kutta", not "Range-Kutta". Two people, one with name "Runge", the other with name "Kutta".

clear all, close all, clc
s = 7.5;                    % semi span
c = 2;                      % chord
hline = 80;                 % hinge line as percentage of chord
xh = hline / 100 * c;       % distance of hinge line aft of leading edge
fa = 40;                    % elastic axis as percentage of chord
xf = fa / 100 * c;          % position of elastic axis aft of le
xac = 0.25 * c;             % distance of ac aft of le
xe = xf - xac;              % distance of elastic axis aft of ac
e = xe / c;                 % eccentricity of fa aft of ac per chord
heave = 3;                  % heave angle
alpha = 3;                  % pitch angle
beta = 3;                   % control surface angle


chat = 2 * (xh / c) - 1;

T4 = - acos(chat) + chat * sqrt(1 - chat ^ 2);
T10 = sqrt(1 - chat ^ 2) + acos(chat);
T11 = acos(chat) * (1 - 2 * chat) + sqrt(1 - chat ^ 2) * (2 - chat);
T12 = sqrt(1 - chat ^ 2) * (2 + chat) - acos(chat) * (2 * chat + 1);

% Aerodynamic coefficients

aw = 2 * pi;                % lift per incidence
ac = aw * T10 / pi;         % lift per control rotation
bw = e * aw;                % pitching moment per control rotation
bc = e * aw * T10 / pi;     % pitching moment per control rotation
cw = - T12 / 2;             % hinge moment per incidence
cc = - T12 * T10 / 2 / pi;  % hinge moment per control rotation

Mtdot = -1.2;               % unsteady torsional damping term
Mbdot = -0.1;               % unsteady control rotation damping term

rho = 1.225;                % air density

% Mass data
m = 400;                    % mass per unit area
mmain = m;                  % mass per unit area of main surface
mcont = m;                  % mass per unit area of control surface

% Stiffness data
EI =  4e+7;                 % flexural rigidity
GJ =  8e+6;                 % torsional rigidity
kbeta = 1e+4;               % control stiffness per unit span
kh = 4 * EI / s^3;          % bending spring stiffness
ka = GJ / s;                % torsion spring stiffness
ka3 = 100*ka*alpha^3;       %non linear stiffness

GJstr = sprintf('%0.5g',GJ); EIstr = sprintf('%0.5g',EI);
astring = ['xf/c = ', num2str(xf/c),' xh/c = ',num2str(xh/c),' EI = ',EIstr,...
    ' GJ = ',GJstr,' kbeta = ', num2str(kbeta)];
disp(astring)

xf/c = 0.4 xh/c = 0.8 EI = 4e+07 GJ = 8e+06 kbeta = 10000

modes = 3;

% Inertia matrix
A = zeros(modes, modes);
% Main surface
Amain = zeros(modes, modes);
Amain(1,1) = xh * s / 5;
Amain(1,2) = (xh^2 / 2 - xh * xf) * s / 4;
Amain(2,1) = Amain(1,2);
Amain(2,2) = (xh^3 / 3 - xh^2 * xf + xf^2 * xh) * s / 3;
% Control surface
Acont = zeros(modes, modes);
Acont(1,1) = (c - xh) * s / 5;
Acont(1,2) = ((c^2 - xh^2) / 2 - xf * (c - xh)) * s / 4;
Acont(1,3) = ((c^2 - xh^2) / 2 - xh * (c - xh)) * s / 3;
Acont(2,2) = ((c^3 - xh^3) / 3 - xf * (c^2 - xh^2) + xf^2 * (c - xh)) * s / 3;
Acont(2,3) = ((c^3 - xh^3) / 3 - xf * (c^2 - xh^2) / 2 - xh * (c^2 - xh^2) / 2 + xf * xh * (c - xh)) * s / 2;
Acont(3,3) = ((c^3 - xh^3) / 3 - xh * (c^2 - xh^2) + xh^2 * (c - xh)) * s;
Acont(2,1) = Acont(1,2);
Acont(3,1) = Acont(1,3);
Acont(3,2) = Acont(2,3);

A = mmain * Amain + mcont * Acont;

% Structural damping matrix is zero
D = zeros(modes, modes);

B = zeros(modes, modes);     % aero damping - based on rho*V*B
C = zeros(modes, modes);     % aero stiffness - based on rho*V^2*C

B(1,1) = aw * c * s / 10;
B(1,2) = 0;
B(1,3) = 0;
B(2,1) = - bw * c^2 * s / 8;
B(2,2) = - Mtdot * c^3 * s / 24;
B(2,3) = 0;
B(3,1) = - cw * c^2 * s / 6;
B(3,2) = 0;
B(3,3) = - Mbdot * c^3 * s / 8;

C(1,1) = 0;
C(1,2) = aw * c * s / 8;
C(1,3) = ac * c * s / 6;
C(2,1) = 0;
C(2,2) = - bw * c^2 * s / 6;
C(2,3) = - bc * c^2 * s / 4;
C(3,1) = 0;
C(3,2) = - cw * c^2 * s / 4;
C(3,3) = - cc * c^2 * s / 2;


E = zeros(modes, modes);
E(1,1) = kh; 
E(2,2) = ka;
E(3,3) = kbeta * s;
n = 3;
V = 117;
t = linspace(0,10,500);
Mat=[zeros(n,n),eye(n,n);
    -A\(rho*V*V*C+E), -A\(rho*V*B+D)];

%X1 = [x(1)
%    x(2)
%    x(3)
%    x(4)
%    x(5)
%    x(6)];

%x(1) = x1;
%x(2) = x2;
%x(3) = x3;
%x(4) = x4;
%x(5) = x5;
%x(6) = x6;

%não linearidade em alpha

%X2 = [0
%    x(2)^3
%    0
%    0
%    0
%    0];

size(Mat)

ans = 1×2

6 6

size([-ka3*A\[0; 1 ;0] zeros(3,1)])

ans = 1×2

3 2

Mat3 = @(t,x) Mat*x + [-ka3*A\[0; 1 ;0] zeros(3,1)]*[0; x(2)^3;0;0;0;0];
x0 = zeros(6,1);
tspan = [0 1];
[T,X] = ode45(Mat3,tspan,x0)

Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To operate on each element of the matrix individually, use TIMES (.*) for elementwise multiplication.

Error in solution>@(t,x)Mat*x+[-ka3*A\[0;1;0],zeros(3,1)]*[0;x(2)^3;0;0;0;0] (line 142)
Mat3 = @(t,x) Mat*x + [-ka3*A\[0; 1 ;0] zeros(3,1)]*[0; x(2)^3;0;0;0;0];

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

plot(T,X)