# Solving Coupled Differential Equations

조회 수: 5(최근 30일)
댓글: 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 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표시숨기기 없음
Thank you, your code helped me figure out what I needed to know about ode45 that I didnt understand before.

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

R2020b

### Community Treasure Hunt

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

Start Hunting!

Translated by