Solution of a 2nd order non linear implicit differential equation using ode15i implicit solver
32 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Parthajit
am 22 Aug. 2024 um 5:10
Kommentiert: Parthajit
am 22 Aug. 2024 um 10:53
This is a nonlinear differential equation of 2nd order and implicit
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1757439/image.png)
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
0 Kommentare
Akzeptierte Antwort
Sam Chak
am 22 Aug. 2024 um 10:30
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
Weitere Antworten (2)
Naga
am 22 Aug. 2024 um 6:49
Bearbeitet: Naga
am 22 Aug. 2024 um 6:55
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
am 22 Aug. 2024 um 10:21
Bearbeitet: Torsten
am 22 Aug. 2024 um 10:21
% 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
0 Kommentare
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!