ODE 45 to solve coupled ODE

1 visualización (últimos 30 días)
aakash dewangan
aakash dewangan el 28 de Mayo de 2021
Comentada: aakash dewangan el 30 de Mayo de 2021
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 comentarios
Jan
Jan el 29 de Mayo de 2021
Which program has created a function including "1.8e1 / 5.0" and (6.0 / 5.0)?!
aakash dewangan
aakash dewangan el 30 de Mayo de 2021
That was the expression of ftotal, when i printed it..

Iniciar sesión para comentar.

Respuestas (1)

Torsten
Torsten el 29 de Mayo de 2021
Editada: Torsten el 29 de Mayo de 2021
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)

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by