Solving Coupled Differential Equations

조회 수: 1(최근 30일)
Braden Kerr
Braden Kerr 2020년 12월 19일
댓글: Braden Kerr 2020년 12월 20일
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.

채택된 답변

Alan Stevens
Alan Stevens 2020년 12월 19일
편집: Alan Stevens 2020년 12월 19일
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
Braden Kerr
Braden Kerr 2020년 12월 20일
Thank you, your code helped me figure out what I needed to know about ode45 that I didnt understand before.

댓글을 달려면 로그인하십시오.

추가 답변(0개)

제품


릴리스

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by