Solving Coupled Differential Equations
3 次查看(过去 30 天)
显示 更早的评论
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.
0 个评论
采纳的回答
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
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!