Solving Coupled Differential Equations

2 views (last 30 days)
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.

Accepted Answer

Alan Stevens
Alan Stevens on 19 Dec 2020
Edited: Alan Stevens on 19 Dec 2020
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
  1 Comment
Braden Kerr
Braden Kerr on 20 Dec 2020
Thank you, your code helped me figure out what I needed to know about ode45 that I didnt understand before.

Sign in to comment.

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!