MATLAB: Solving Coupled Differential Equations

coupled equationsMATLABode45

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 4×1 matrix and f(x) is also a 4×1. 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.

Best Answer

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
Related Question