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

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

Answer 1

[Sam Chak]

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

[Parthajit]

Thank you so much Sir. It helped me a lot.

Answer 2

[Naga]

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))

[Parthajit]

Can you please provide the whole code which peovides the velocity history and time history

Answer 3

[Torsten]

% 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