ODE 45 to solve coupled ODE
2 次查看(过去 30 天)
显示 更早的评论
Hello,
I have written a code to solve a system of coupled second order differential equations using ODE 45.
I wrote this code by modifying the available code. I want to know that what are the initial conditions being applied to which variable (p1, p2, p3). Please let me know what are the initial conditions for p1, Dp1, p2, Dp2, p3, Dp3. And also which column of "y" in "pSol" corresponds to p1, p2 and p3 (in other words how can i get p1, p2, p3 from pSol).
Thanks in advance,
clc; clear all; close all
syms p1(t) p2(t) p3(t)
rho L m v T k G
% Parameters
rho = 9000; T = 15000; L = 100; m = 5; v = 120; k = 2000; G = 0.2
Dp1 = diff(p1); D2p1 = diff(p1,2); Dp2 = diff(p2); D2p2 = diff(p2,2); Dp3 = diff(p3); D2p3 = diff(p3,2);
% Mass matrix terms
AA = rho*L/2 + m*(sin(pi*v*t/L))^2;
BB = m*sin(2*pi*v*t/L)*sin(pi*v*t/L);
CC = m*sin(3*pi*v*t/L)*sin(pi*v*t/L);
DD = rho*L/2 + m*(sin(2*pi*v*t/L))^2;
EE = m*sin(2*pi*v*t/L)*sin(3*pi*v*t/L);
FF = rho*L/2 + m*(sin(3*pi*v*t/L))^2;
% Stiffness matrix terms
GG = T*(pi/L)^2*(L/2) + k*(sin(pi*v*t/L))^2;
HH = k*sin(2*pi*v*t/L)*sin(pi*v*t/L);
II = k*sin(pi*v*t/L)*sin(3*pi*v*t/L);
JJ = T*(2*pi/L)^2*(L/2) + k*(sin(2*pi*v*t/L))^2;
KK = k*sin(2*pi*v*t/L)*sin(3*pi*v*t/L);
LL = T*(3*pi/L)^2*(L/2) + k*(sin(3*pi*v*t/L))^2;
% RHS
MM = k*G*sin(pi*v*t/L);
NN = k*G*sin(2*pi*v*t/L);
OO = k*G*sin(3*pi*v*t/L);
Eq1 = AA*diff(p1,t,2) + BB*diff(p2,t,2) + CC*diff(p3,t,2) + GG*p1 + HH*p2 + II*p3 == MM;
Eq2 = BB*diff(p1,t,2) + DD*diff(p2,t,2) + EE*diff(p3,t,2) + HH*p1 + JJ*p2 + KK*p3 == NN;
Eq3 = CC*diff(p1,t,2) + EE*diff(p2,t,2) + FF*diff(p3,t,2) + II*p1 + KK*p2 + LL*p3 == OO;
[V,S] = odeToVectorField(Eq1, Eq2, Eq3);
ftotal = matlabFunction(V, 'Vars',{'t','Y'});
interval = [0 L/v];
y0 = [1 2; 3 4; 5 6]; %initial conditions
% v-k2
pSol = ode45( @(t,Y)ftotal(t,Y),interval,y0);
4 个评论
回答(1 个)
Torsten
2021-5-29
编辑:Torsten
2021-5-29
To be honest, I'd prefer to know what Matlab solves.
Order the variables as
[z1,z2,z3,z4,z5,z6] = [y1,y2,y3,y1',y2',y3']
Then your adapted mass matrix becomes
function MM = mass(t,z)
%Define M
MM = [eye(3),zeros(3);zeros(3),M];
end
and your right-hand side vector becomes
function dz = fun(t,z)
% Define K and F
dz = [zeros(3),eye(3);-K,zeros(3)]*z + [F;zeros(3,1)]
end
and the call to ode15s
options = odeset('Mass',@mass);
[T,Z] = ode15s(@fun,tspan,z0,options)
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!