Ordinary Differential Equations¶
Introduction¶
In [1]:
f = @(t, x) -x.^2 + t;
In [2]:
tf = 5.0;
N = 30;
t = linspace(-2, tf, N);
x = linspace(-4, 4, N);
[T, X] = meshgrid(t, x);
G = f(T, X); % G for gradient
ANG = atan(G);
U = ones(size(T)) .* cos(ANG);
V = ones(size(T)) .* sin(ANG);
figure(1, 'Position', [0, 0, 1500, 1000]);
quiver(T - U / 2, X - V / 2, U, V);
hold on;
% draw x-axis
plot(0:tf, zeros(size(0:tf)), 'Color', 'black', 'LineWidth', 5);
% draw y-axis
plot(zeros(size(0:tf)), 0:tf, 'Color', 'black', 'LineWidth', 5);
xlim([-2.1, tf + 0.1]); % actually t - time
xlabel('time, t', 'FontSize', 25);
ylim([-4.1, 4.1]); % actually x - trajectory
ylabel('trajectory, x', 'FontSize', 25);
In [3]:
tf = 5.0;
figure(2, 'Position', [0, 0, 1500, 1000]);
quiver(T - U / 2, X - V / 2, U, V);
hold on;
% draw x-axis
plot(0:tf, zeros(size(0:tf)), 'Color', 'black', 'LineWidth', 5);
% draw y-axis
plot(zeros(size(0:tf)), 0:tf, 'Color', 'black', 'LineWidth', 5);
xlim([-2.1, tf + 0.1]); % actually t - time
xlabel('time, t', 'FontSize', 25);
ylim([-4.1, 4.1]); % actually x - trajectory
ylabel('trajectory, x', 'FontSize', 25);
% ---------------------------------------------------
x0 = 0.0;
t = linspace(0, tf, 1e3);
x = zeros(1, length(t));
x(1) = x0;
dt = (t(end) - t(1)) / (length(t) - 1);
for i = 1:(length(t) - 1)
x(i + 1) = x(i) + dt * f(t(i), x(i));
end
plot(t, x, 'LineWidth', 10, 'DisplayName', 'x0 = 0.0');
% ---------------------------------------------------
x0 = 3.0;
t = linspace(0, tf, 1e3);
x = zeros(1, length(t));
x(1) = x0;
dt = (t(end) - t(1)) / (length(t) - 1);
for i = 1:(length(t) - 1)
x(i + 1) = x(i) + dt * f(t(i), x(i));
end
plot(t, x, 'LineWidth', 10, 'DisplayName', 'x0 = 3.0');
% ---------------------------------------------------
x0 = -0.5;
t = linspace(0, tf, 1e3);
x = zeros(1, length(t));
x(1) = x0;
dt = (t(end) - t(1)) / (length(t) - 1);
for i = 1:(length(t) - 1)
x(i + 1) = x(i) + dt * f(t(i), x(i));
end
plot(t, x, 'LineWidth', 10, 'DisplayName', 'x0 = -0.5');
% ---------------------------------------------------
x0 = -0.8;
t = linspace(0, tf, 1e3);
x = zeros(1, length(t));
x(1) = x0;
dt = (t(end) - t(1)) / (length(t) - 1);
for i = 1:(length(t) - 1)
x(i + 1) = x(i) + dt * f(t(i), x(i));
end
plot(t, x, 'LineWidth', 10, 'DisplayName', 'x0 = -0.7');
lgd = legend('show');
set(lgd, 'FontSize', 20);
In [ ]:
Higher Order¶
Harmonic Oscillator with damping
$\ddot{x} = -0.3 \dot{x} - x \implies \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = \begin{bmatrix} x_2 \\ -0.3 x_2 - x_1 \end{bmatrix} $
In [1]:
x0 = [5.0, 0.0];
f = @(t, x) [x(2), -0.5 * x(2) - x(1)];
t = linspace(0, 10, 1e3);
x = zeros(length(t), 2);
x(1, :) = x0;
dt = (t(end) - t(1)) / (length(t) - 1);
for i = 1:(length(t) - 1)
x(i + 1, :) = x(i, :) + dt * f(t, x(i, :));
end
figure(3);
plot(t, x(:, 1), 'LineWidth', 10, 'DisplayName', 'position, x_1');
hold on;
plot(t, x(:, 2), 'LineWidth', 10, 'DisplayName', 'speed, x_2');
xlabel('Time, t', 'FontSize', 18);
ylabel('State, x', 'FontSize', 18);
legend('show');
Using in-built solvers¶
Harmonic Oscillator with damping
$\ddot{x} = -0.3 \dot{x} - x \implies \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = \begin{bmatrix} x_2 \\ -0.3 x_2 - x_1 \end{bmatrix} $
In [2]:
x0 = [5.0, 0.0];
f = @(t, x) [x(2), -0.5 * x(2) - x(1)];
t = linspace(0, 10, 1e3);
[t, x] = ode45(@(t, x) f(t, x), t, x0);
figure(3);
plot(t, x(:, 1), 'LineWidth', 10, 'DisplayName', 'position, x_1');
hold on;
plot(t, x(:, 2), 'LineWidth', 10, 'DisplayName', 'speed, x_2');
xlabel('Time, t', 'FontSize', 18);
ylabel('State, x', 'FontSize', 18);
Mass on a spring¶
$m a = -k x$
$m \ddot{x} = -k x$
$\ddot{x} = -\frac{k}{m} x$
$\begin{bmatrix} \dot{x_1} \\ \dot{x_2} \\ \end{bmatrix} = \begin{bmatrix} x_2 \\ -\frac{k}{m} x_1 \\ \end{bmatrix}$
In [3]:
function du = f(t, u, k, m)
x = u(1);
v = u(2);
a = -k / m * x;
du = [v; a];
end
In [5]:
u0 = [5, 0];
tspan = linspace(0, 10, 1e2);
k = 3;
m = 1;
[T, U] = ode45(@(t, u) f(t, u, k, m), tspan, u0);
figure();
% position, x
plot(T, U(:, 1), 'LineWidth', 7, 'DisplayName', 'x');
hold on;
% velocity, v
plot(T, U(:, 2), 'LineWidth', 7, 'DisplayName', 'v');
xlabel('Time, t', 'FontSize', 18);
ylabel('State, u_i', 'FontSize', 18);
set(gca(), 'FontSize', 18);
lgd = legend('show');
set(lgd, 'FontSize', 18);
In [ ]: