Solution of a 2nd order non linear implicit differential equation using ode15i implicit solver

8 次查看（过去 30 天）
Parthajit
Parthajit 2024-8-22
评论： Parthajit 2024-9-9
This is a nonlinear differential equation of 2nd order and implicit
I am providing the code I have tried but Iam not getting hoe the velocity changes with time at different time what is the value of velocity.
% Initial conditions
y0 = [0 2]; % initial displacement and velocity
yp0_guess = [0; 0]; % Initial guess for derivatives
% Time span
tspan = [0 20]; % time span for the solution
% Solve using ode15i
[t, y] = ode15i(@implicitODE, tspan, y0, yp0_guess);
% Plot the displacement over time
figure;
subplot(2,1,1);
plot(t, y(:,1));
xlabel('Time (s)');
ylabel('Displacement (y)');
title('Displacement vs Time');
% Plot the velocity over time
subplot(2,1,2);
plot(t,y(:,2)); % y(:,2) corresponds to the velocity (dy/dt)
xlabel('Time (s)');
ylabel('Velocity (dy/dt)');
title('Velocity vs Time');
function res = implicitODE(t, y, yp)
% Constants
pi = 3.141592653589793;
k1 = 4.4049e-18;
k2 = 0.1101;
a = 17.71e-3;
% Extract variables
y1 = y(1); % y1 corresponds to y(t)
y2 = y(2); % y2 corresponds to dy/dt
yp1 = yp(1); % yp1 corresponds to y2
yp2 = yp(2); % yp2 corresponds to y2'
% Compute the residuals
res = zeros(2,1);
res(1) = yp1 - y2; % y1' = y2
res(2) = yp2 - (-k1 * y1^2 * y2 / (a^2 + y1^2)^2.5 / k2)-9.81; % y2' equation
end

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

采纳的回答

Sam Chak
Sam Chak 2024-8-22
Ran in:
Hi @Parthajit
Please check if the following is what you expected from the simulation.
tspan = [0 20];
y0 = [0; 2];
[t, y] = ode45(@ode, tspan, y0);
plot(t, y), grid on
xlabel('t')
legend('displacement', 'velocity', 'location', 'south')
function dydt = ode(t, y)
% Constants
k1 = 4.4049e-18;
k2 = 0.1101;
a = 17.71e-3;
% System
dydt = zeros(2, 1);
dydt(1) = y(2);
dydt(2) = - (k1*(y(1)^2)*y(2))/(k2*(a^2 + y(1)^2)^(5/2)) - 9.81;
end

更多回答（2 个）

Naga
Naga 2024-8-22
编辑：Naga 2024-8-22
Hi Parthajith,
I understand you want to know the velocity values at different times. To find that, you can read the values of T and y(:,2), where T represents the time stamps and y(:,2) represents the velocity at those particular time stamps. I'm attaching the values below:
T = [0, 0.0000, 0.0001, 0.0001, 0.0001, 0.0002, 0.0002, 0.0003, 0.0004, 0.0005, ...
0.0007, 0.0011, 0.0018, 0.0034, 0.0066, 0.0130, 0.0257, 0.0512, 0.1022, 0.2041, ...
0.4080, 0.8157, 1.6311, 3.2620, 5.2620, 7.2620, 9.2620, 11.2620, 13.2620, 15.2620, ...
17.2620, 19.2620, 20.0000]; %The array T contains the time stamps at which the velocities are recorded.
y(:,2) = [2.0000, 2.0003, 2.0005, 2.0008, 2.0010, 2.0015, 2.0020, 2.0025, 2.0035, 2.0044, ...
2.0064, 2.0103, 2.0181, 2.0337, 2.0650, 2.1275, 2.2525, 2.5025, 3.0024, 4.0023, ...
6.0022, 10.0018, 18.0011, 33.9998, 53.6198, 73.2398, 92.8598, 112.4798, 132.0998, ...
151.7198, 171.3398, 190.9598, 198.2000]; %The array y(:,2) represents the velocity values corresponding to each time stamp in T
plot(T,y(:,2))
Torsten
Torsten 2024-8-22
编辑：Torsten 2024-8-22
Ran in:
% Initial conditions
y0 = [0 2]; % initial displacement and velocity
%yp0_guess = [0; 0]; % Initial guess for derivatives
% Time span
tspan = [0 20]; % time span for the solution
% Solve using ode15i
%[t, y] = ode15i(@implicitODE, tspan, y0, yp0_guess);
[t, y] = ode15s(@nonimplicitODE, tspan, y0);
% Plot the displacement over time
figure;
subplot(2,1,1);
plot(t, y(:,1));
xlabel('Time (s)');
ylabel('Displacement (y)');
title('Displacement vs Time');
% Plot the velocity over time
subplot(2,1,2);
plot(t,y(:,2)); % y(:,2) corresponds to the velocity (dy/dt)
xlabel('Time (s)');
ylabel('Velocity (dy/dt)');
title('Velocity vs Time');
%function res = implicitODE(t, y, yp)
function dy = nonimplicitODE(t, y)
% Constants
pi = 3.141592653589793;
k1 = 4.4049e-18;
k2 = 0.1101;
a = 17.71e-3;
% Extract variables
y1 = y(1); % y1 corresponds to y(t)
y2 = y(2); % y2 corresponds to dy/dt
%yp1 = yp(1); % yp1 corresponds to y2
%yp2 = yp(2); % yp2 corresponds to y2'
% Compute the residuals
%res = zeros(2,1);
%res(1) = yp1 - y2; % y1' = y2
%res(2) = yp2 - (-k1 * y1^2 * y2 / (a^2 + y1^2)^2.5 / k2)-9.81; % y2' equation
dy = zeros(2,1);
dy(1) = y2;
dy(2) = -k1 * y1^2 * y2 / ( k2*(a^2 + y1^2)^2.5 ) - 9.81;
end

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by