Solving Coupled Differential Equations

Hello,

I am trying to model the flutter in a wing and have two coupled equations of motion. There are two parameters I am modeling, h and alpha. Equation 1 is a function of h, h', h'', alpha' and alpha''. The second is a function of alpha'', h'', alpha', and a. Through linearization I can reduce such a system to something of the form x' = f(x) where x' is a 4x1 matrix and f(x) is also a 4x1. I am attempting to solve such a condition with given values of h, alpha, h' and alpha', but I am unsure how to input these into ode45 properly. The full code I currently have is as follows.

% Given parameters
syms h a h_dot a_dot
m = 1;
x_m = 0.05;
e = 0.4;
zeta = 0.1;
w_h = 0.5;
w_a = 1.0;
I_a = 0.25;
D_a = 0.2;


%GUESS U 
U = 5;
%

L = D_a*U^2*(a+h_dot/U);
M = [1 x_m; (m*x_m)/I_a 1];
C = [2*zeta*w_h 0; 0 2*zeta*w_a];
K = [w_h^2 0; 0 w_a^2];
G1 = [-L/m; (L*e)/I_a];
G2 = [0;0];
I = [1 0; 0 1];
Z = [0 0; 0 0];
A = [0.5*C M; I Z];
B = [K 0.5*C; Z -I];
F = [G1; G2];

y = [h; a];
y_dot = [h_dot; a_dot];
x_vector = [y; y_dot];

Ainv = inv(A);

%RHS = -Ainv*B*x + Ainv*F;
% Initial Conditions
h = 0;
a = 0;
h_dot = 0;
a_dot = 0;

time = [0 5];

fun = @(time,x) [x(2); -Ainv*B*x(1) + Ainv*F];

[T,X] = ode45(fun, time, x_vector);

plot(T, X(:,1));
hold on 
plot(T, X(:,2));
hold off

I provide this to answer any questions about variables or other names. However, the important section is below

%RHS = -Ainv*B*x + Ainv*F;
% Initial Conditions
h = 0;
a = 0;
h_dot = 0;
a_dot = 0;

time = [0 5];

fun = @(time,x) [x(2); -Ainv*B*x(1) + Ainv*F];

[T,X] = ode45(fun, time, x_vector);

plot(T, X(:,1));
hold on 
plot(T, X(:,2));
hold off

I get quite a few errors and unsure how to resolve them, any adivce is appreciated, thank you.

Answer 1

[Alan Stevens]

More like this perhaps, though the end result is rather boring as your forcing function, F, is all zeros and your initial conditions are all zero!

% A*x' + B*x = F
% x' = A\(-B*x+F)

% Initial Conditions
h = 0;
a = 0;
h_dot = 0;
a_dot = 0;
 
x_vector = [h; a; h_dot; a_dot];
time = [0 5];

[T,X] = ode45(@fun, time, x_vector);

plot(T, X(:,1));
hold on 
plot(T, X(:,2));
hold off

function dXdt = fun(~, x)
        m = 1;
        x_m = 0.05;
        e = 0.4;
        zeta = 0.1;
        w_h = 0.5;
        w_a = 1.0;
        I_a = 0.25;
        D_a = 0.2;
        U = 5;
        
        a = x(2); h_dot = x(3);
        
        L = D_a*U^2*(a+h_dot/U);
        M = [1 x_m; (m*x_m)/I_a 1];
        C = [2*zeta*w_h 0; 0 2*zeta*w_a];
        K = [w_h^2 0; 0 w_a^2];
        G1 = [-L/m; (L*e)/I_a];
        G2 = [0;0];
        I = [1 0; 0 1];
        Z = [0 0; 0 0];
        A = [0.5*C M; I Z];
        B = [K 0.5*C; Z -I];
        F = [G1; G2];

        dXdt = A\(-B*x + F);

end

[Braden Kerr]

Thank you, your code helped me figure out what I needed to know about ode45 that I didnt understand before.