ODE 45 to solve coupled ODE

2021 5 28
1 个回答
更新时间：2021 5 30
8 次查看（30 天）

0 个投票

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 个评论
显示 2更早的评论 隐藏 2更早的评论

Torsten
Torsten 2021-5-28
Show us ftotal.
aakash dewangan
aakash dewangan 2021-5-29
Torsten, ftotal expression in command window looks like this
@(t,Y)[Y(2);(sin(pi.*t.*(1.2e1./5.0)).*8.0e1-sin(pi.*t.*(6.0./5.0)).^2.*Y(1).*6.579736267392906e-3-sin(pi.*t.*(1.2e1./5.0)).^2.*Y(1).*4.0e2-sin(pi.*t.*(1.8e1./5.0)).^2.*Y(1).*6.579736267392906e-3-Y(1).*5.921762640653616e2-sin(pi.*t.*(6.0./5.0)).*sin(pi.*t.*(1.2e1./5.0)).*Y(3).*3.999983550659332e2-sin(pi.*t.*(1.2e1./5.0)).*sin(pi.*t.*(1.8e1./5.0)).*Y(5).*3.999851955933984e2)./(sin(pi.*t.*(6.0./5.0)).^2+sin(pi.*t.*(1.2e1./5.0)).^2+sin(pi.*t.*(1.8e1./5.0)).^2+9.0e4);Y(4);(sin(pi.*t.*(6.0./5.0)).*8.0e1-sin(pi.*t.*(6.0./5.0)).^2.*Y(3).*4.0e2-sin(pi.*t.*(1.2e1./5.0)).^2.*Y(3).*1.644934066848227e-3-sin(pi.*t.*(1.8e1./5.0)).^2.*Y(3).*1.644934066848227e-3-Y(3).*1.480440660163404e2-sin(pi.*t.*(6.0./5.0)).*sin(pi.*t.*(1.2e1./5.0)).*Y(1).*3.999934202637326e2-sin(pi.*t.*(6.0./5.0)).*sin(pi.*t.*(1.8e1./5.0)).*Y(5).*3.999851955933984e2)./(sin(pi.*t.*(6.0./5.0)).^2+sin(pi.*t.*(1.2e1./5.0)).^2+sin(pi.*t.*(1.8e1./5.0)).^2+9.0e4);Y(6);(sin(pi.*t.*(1.8e1./5.0)).*8.0e1-sin(pi.*t.*(6.0./5.0)).^2.*Y(5).*1.480440660163404e-2-sin(pi.*t.*(1.2e1./5.0)).^2.*Y(5).*1.480440660163404e-2-sin(pi.*t.*(1.8e1./5.0)).^2.*Y(5).*4.0e2-Y(5).*1.332396594147063e3-sin(pi.*t.*(6.0./5.0)).*sin(pi.*t.*(1.8e1./5.0)).*Y(3).*3.999983550659332e2-sin(pi.*t.*(1.2e1./5.0)).*sin(pi.*t.*(1.8e1./5.0)).*Y(1).*3.999934202637326e2)./(sin(pi.*t.*(6.0./5.0)).^2+sin(pi.*t.*(1.2e1./5.0)).^2+sin(pi.*t.*(1.8e1./5.0)).^2+9.0e4)]
Jan
Jan 2021-5-29
Which program has created a function including "1.8e1 / 5.0" and (6.0 / 5.0)?!
aakash dewangan
aakash dewangan 2021-5-30
That was the expression of ftotal, when i printed it..

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2021-5-29
编辑：Torsten 2021-5-29

0 个投票

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)

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by