Main Content
Results for
Following on from my previous post The Non-Chaotic Duffing Equation, now we will study the chaotic behaviour of the Duffing Equation
P.s:Any comments/advice on improving the code is welcome.
The Original Duffing Equation is the following:
Let 
. This implies that
. This implies thatThen we rewrite it as a System of First-Order Equations
Using the substitution 
for 
, the second-order equation can be transformed into the following system of first-order equations:
for , the second-order equation can be transformed into the following system of first-order equations:
Exploring the Effect of γ.
% Define parameters
gamma = 0.1;
alpha = -1;
beta = 1;
delta = 0.1;
omega = 1.4;
% Define the system of equations
odeSystem = @(t, y) [y(2);
-delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];
% Initial conditions
y0 = [0; 0]; % x(0) = 0, v(0) = 0
% Time span
tspan = [0 200];
% Solve the system
[t, y] = ode45(odeSystem, tspan, y0);
% Plot the results
figure;
plot(t, y(:, 1));
xlabel('Time');
ylabel('x(t)');
title('Solution of the nonlinear system');
grid on;
% Plot the phase portrait
figure;
plot(y(:, 1), y(:, 2));
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait');
grid on;
% Define the tail (e.g., last 10% of the time interval)
tail_start = floor(0.9 * length(t)); % Starting index for the tail
tail_end = length(t); % Ending index for the tail
% Plot the tail of the solution
figure;
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), 'r', 'LineWidth', 1.5);
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait - Tail of the Solution');
grid on;
% Define parameters
gamma = 0.318;
alpha = -1;
beta = 1;
delta = 0.1;
omega = 1.4;
% Define the system of equations
odeSystem = @(t, y) [y(2);
-delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];
% Initial conditions
y0 = [0; 0]; % x(0) = 0, v(0) = 0
% Time span
tspan = [0 800];
% Solve the system
[t, y] = ode45(odeSystem, tspan, y0);
% Plot the results
figure;
plot(t, y(:, 1));
xlabel('Time');
ylabel('x(t)');
title('Solution of the nonlinear system');
grid on;
% Plot the phase portrait
figure;
plot(y(:, 1), y(:, 2));
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait');
grid on;
% Define the tail (e.g., last 10% of the time interval)
tail_start = floor(0.9 * length(t)); % Starting index for the tail
tail_end = length(t); % Ending index for the tail
% Plot the tail of the solution
figure;
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), 'r', 'LineWidth', 1.5);
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait - Tail of the Solution');
grid on;
% Define parameters
gamma = 0.338;
alpha = -1;
beta = 1;
delta = 0.1;
omega = 1.4;
% Define the system of equations
odeSystem = @(t, y) [y(2);
-delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];
% Initial conditions
y0 = [0; 0]; % x(0) = 0, v(0) = 0
% Time span with more points for better resolution
tspan = linspace(0, 200,2000); % Increase the number of points
% Solve the system
[t, y] = ode45(odeSystem, tspan, y0);
% Plot the results
figure;
plot(t, y(:, 1));
xlabel('Time');
ylabel('x(t)');
title('Solution of the nonlinear system');
grid on;
% Plot the phase portrait
figure;
plot(y(:, 1), y(:, 2));
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait');
grid on;
% Define the tail (e.g., last 10% of the time interval)
tail_start = floor(0.9 * length(t)); % Starting index for the tail
tail_end = length(t); % Ending index for the tail
% Plot the tail of the solution
figure;
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), 'r', 'LineWidth', 1.5);
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait - Tail of the Solution');
grid on;
ax = gca;
chart = ax.Children(1);
datatip(chart,0.5581,-0.1126);
% Define parameters
gamma = 0.35;
alpha = -1;
beta = 1;
delta = 0.1;
omega = 1.4;
% Define the system of equations
odeSystem = @(t, y) [y(2);
-delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];
% Initial conditions
y0 = [0; 0]; % x(0) = 0, v(0) = 0
% Time span with more points for better resolution
tspan = linspace(0, 400,3000); % Increase the number of points
% Solve the system
[t, y] = ode45(odeSystem, tspan, y0);
% Plot the results
figure;
plot(t, y(:, 1));
xlabel('Time');
ylabel('x(t)');
title('Solution of the nonlinear system');
grid on;
% Plot the phase portrait
figure;
plot(y(:, 1), y(:, 2));
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait');
grid on;
% Define the tail (e.g., last 10% of the time interval)
tail_start = floor(0.9 * length(t)); % Starting index for the tail
tail_end = length(t); % Ending index for the tail
% Plot the tail of the solution
figure;
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), 'r', 'LineWidth', 1.5);
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait - Tail of the Solution');
grid on;
Studying the attached document Duffing Equation from the University of Colorado, I noticed that there is an analysis of The Non-Chaotic Duffing Equation and all the graphs were created with Matlab. And since the code is not given I took the initiative to try to create the same graphs with the following code.
- Plotting the Potential Energy and Identifying Extrema
% Define the range of x values
x = linspace(-2, 2, 1000);
% Define the potential function V(x)
V = -x.^2 / 2 + x.^4 / 4;
% Plot the potential function
figure;
plot(x, V, 'LineWidth', 2);
hold on;
% Mark the minima at x = ±1
plot([-1, 1], [-1/4, -1/4], 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'g');
% Add LaTeX title and labels
title('Duffing Potential Energy: $$V(x) = -\frac{x^2}{2} + \frac{x^4}{4}$$', 'Interpreter', 'latex');
xlabel('$$x$$', 'Interpreter', 'latex');
ylabel('$$V(x)$$','Interpreter', 'latex');
grid on;
hold off;
- Solving and Plotting the Duffing Equation
% Define the system of ODEs for the non-chaotic Duffing equation
duffing_ode = @(t, X) [X(2);
X(1) - X(1).^3];
% Time span for the simulation
tspan = [0 10];
% Initial conditions [x(0), v(0)]
initial_conditions = [1; 1];
% Solve the ODE using ode45
[t, X] = ode45(duffing_ode, tspan, initial_conditions);
% Extract displacement (x) and velocity (v)
x = X(:, 1);
v = X(:, 2);
% Plot both x(t) and v(t) in the same figure
figure;
plot(t, x, 'b-', 'LineWidth', 2); % Plot x(t) with blue line
hold on;
plot(t, v, 'r--', 'LineWidth', 2); % Plot v(t) with red dashed line
% Add title, labels, and legend
title(' Component curve solutions to $$\ddot{x}-x+x^3=0$$','Interpreter', 'latex');
xlabel('t','Interpreter', 'latex');
ylabel('$$x(t) $$ and $$v(t) $$','Interpreter', 'latex');
legend('$$x(t)$$', ' $$v(t)$$','Interpreter', 'latex');
grid on;
hold off;
% Phase portrait with nullclines, equilibria, and vector field
figure;
hold on;
% Plot phase portrait
plot(x, v,'r', 'LineWidth', 2);
% Plot equilibrium points
plot([0 1 -1], [0 0 0], 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'g');
% Create a grid of points for the vector field
[x_vals, v_vals] = meshgrid(linspace(-2, 2, 20), linspace(-1, 1, 20));
% Compute the vector field components
dxdt = v_vals;
dvdt = x_vals - x_vals.^3;
% Plot the vector field
quiver(x_vals, v_vals, dxdt, dvdt, 'b');
% Set axis limits to [-1, 1]
xlim([-1.7 1.7]);
ylim([-1 1]);
% Labels and title
title('Phase-Plane solutions to $$\ddot{x}-x+x^3=0$$','Interpreter', 'latex');
xlabel('$$ (x)$$','Interpreter', 'latex');
ylabel('$$v(v)$$','Interpreter', 'latex');
grid on;
hold off;
An attractor is called strange if it has a fractal structure, that is if it has non-integer Hausdorff dimension. This is often the case when the dynamics on it are chaotic, but strange nonchaotic attractors also exist. If a strange attractor is chaotic, exhibiting sensitive dependence on initial conditions, then any two arbitrarily close alternative initial points on the attractor, after any of various numbers of iterations, will lead to points that are arbitrarily far apart (subject to the confines of the attractor), and after any of various other numbers of iterations will lead to points that are arbitrarily close together. Thus a dynamic system with a chaotic attractor is locally unstable yet globally stable: once some sequences have entered the attractor, nearby points diverge from one another but never depart from the attractor.
The term strange attractor was coined by David Ruelle and Floris Takens to describe the attractor resulting from a series of bifurcations of a system describing fluid flow. Strange attractors are often differentiable in a few directions, but some are like a Cantor dust, and therefore not differentiable. Strange attractors may also be found in the presence of noise, where they may be shown to support invariant random probability measures of Sinai–Ruelle–Bowen type.
Lorenz
% Lorenz Attractor Parameters
sigma = 10;
beta = 8/3;
rho = 28;
% Lorenz system of differential equations
f = @(t, a) [-sigma*a(1) + sigma*a(2);
rho*a(1) - a(2) - a(1)*a(3);
-beta*a(3) + a(1)*a(2)];
% Time span
tspan = [0 100];
% Initial conditions
a0 = [1 1 1];
% Solve the system using ode45
[t, a] = ode45(f, tspan, a0);
% Plot using scatter3 with time-based color mapping
figure;
scatter3(a(:,1), a(:,2), a(:,3), 5, t, 'filled'); % 5 is the marker size
title('Lorenz Attractor');
xlabel('x(t)');
ylabel('y(t)');
zlabel('z(t)');
grid on;
colorbar; % Add a colorbar to indicate the time mapping
view(3); % Set the view to 3D
Sprott
% Define the parameters
a = 2.07;
b = 1.79;
% Define the system of differential equations
dynamics = @(t, X) [ ...
X(2) + a * X(1) * X(2) + X(1) * X(3); % dx/dt
1 - b * X(1)^2 + X(2) * X(3); % dy/dt
X(1) - X(1)^2 - X(2)^2 % dz/dt
];
% Initial conditions
X0 = [0.63; 0.47; -0.54];
% Time span
tspan = [0 100];
% Solve the system using ode45
[t, X] = ode45(dynamics, tspan, X0);
% Plot the results with color gradient
figure;
colormap(jet); % Set the colormap
c = linspace(1, 10, length(t)); % Color data based on time
% Create a 3D line plot with color based on time
for i = 1:length(t)-1
plot3(X(i:i+1,1), X(i:i+1,2), X(i:i+1,3), 'Color', [0 0.5 0.9]*c(i)/10, 'LineWidth', 1.5);
hold on;
end
% Set plot properties
title('Sprott Attractor');
xlabel('x(t)');
ylabel('y(t)');
zlabel('z(t)');
grid on;
colorbar; % Add a colorbar to indicate the time mapping
view(3); % Set the view to 3D
hold off;
Rössler
% Define the parameters
a = 0.2;
b = 0.2;
c = 5.7;
% Define the system of differential equations
dynamics = @(t, X) [ ...
-(X(2) + X(3)); % dx/dt
X(1) + a * X(2); % dy/dt
b + X(3) * (X(1) - c) % dz/dt
];
% Initial conditions
X0 = [10.0; 0.00; 10.0];
% Time span
tspan = [0 100];
% Solve the system using ode45
[t, X] = ode45(dynamics, tspan, X0);
% Plot the results
figure;
scatter3(X(:,1), X(:,2), X(:,3), 5, t, 'filled');
title('Rössler Attractor');
xlabel('x(t)');
ylabel('y(t)');
zlabel('z(t)');
grid on;
colorbar; % Add a colorbar to indicate the time mapping
view(3); % Set the view to 3D
Rabinovich-Fabrikant
%% Parameters for Rabinovich-Fabrikant Attractor
alpha = 0.14;
gamma = 0.10;
dt = 0.01;
num_steps = 5000;
% Initial conditions
x0 = -1;
y0 = 0;
z0 = 0.5;
% Preallocate arrays for performance
x = zeros(1, num_steps);
y = zeros(1, num_steps);
z = zeros(1, num_steps);
% Set initial values
x(1) = x0;
y(1) = y0;
z(1) = z0;
% Generate the attractor
for i = 1:num_steps-1
x(i+1) = x(i) + dt * (y(i)*(z(i) - 1 + x(i)^2) + gamma*x(i));
y(i+1) = y(i) + dt * (x(i)*(3*z(i) + 1 - x(i)^2) + gamma*y(i));
z(i+1) = z(i) + dt * (-2*z(i)*(alpha + x(i)*y(i)));
end
% Create a time vector for color mapping
t = linspace(0, 100, num_steps);
% Plot using scatter3
figure;
scatter3(x, y, z, 5, t, 'filled'); % 5 is the marker size
title('Rabinovich-Fabrikant Attractor');
xlabel('x(t)');
ylabel('y(t)');
zlabel('z(t)');
grid on;
colorbar; % Add a colorbar to indicate the time mapping
view(3); % Set the view to 3D
References
This project discusses predator-prey system, particularly the Lotka-Volterra equations,which model the interaction between two sprecies: prey and predators. Let's solve the Lotka-Volterra equations numerically and visualize the results.% Define parameters
% Define parameters
alpha = 1.0; % Prey birth rate
beta = 0.1; % Predator success rate
gamma = 1.5; % Predator death rate
delta = 0.075; % Predator reproduction rate
% Define the symbolic variables
syms R W
% Define the equations
eq1 = alpha * R - beta * R * W == 0;
eq2 = delta * R * W - gamma * W == 0;
% Solve the equations
equilibriumPoints = solve([eq1, eq2], [R, W]);
% Extract the equilibrium point values
Req = double(equilibriumPoints.R);
Weq = double(equilibriumPoints.W);
% Display the equilibrium points
equilibriumPointsValues = [Req, Weq]
% Solve the differential equations using ode45
lotkaVolterra = @(t,Y)[alpha*Y(1)-beta*Y(1)*Y(2);
delta*Y(1)*Y(2)-gamma*Y(2)];
% Initial conditions
R0 = 40;
W0 = 9;
Y0 = [R0, W0];
tspan = [0, 100];
% Solve the differential equations
[t, Y] = ode45(lotkaVolterra, tspan, Y0);
% Extract the populations
R = Y(:, 1);
W = Y(:, 2);
% Plot the results
figure;
subplot(2,1,1);
plot(t, R, 'r', 'LineWidth', 1.5);
hold on;
plot(t, W, 'b', 'LineWidth', 1.5);
xlabel('Time (months)');
ylabel('Population');
legend('R', 'W');
grid on;
subplot(2,1,2);
plot(R, W, 'k', 'LineWidth', 1.5);
xlabel('R');
ylabel('W');
grid on;
hold on;
plot(Req, Weq, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');
legend('Phase Trajectory', 'Equilibrium Point');
Now, we need to handle a modified version of the Lotka-Volterra equations. These modified equations incorporate logistic growth fot the prey population.
These equations are:
% Define parameters
alpha = 1.0;
K = 100; % Carrying Capacity of the prey population
beta = 0.1;
gamma = 1.5;
delta = 0.075;
% Define the symbolic variables
syms R W
% Define the equations
eq1 = alpha*R*(1 - R/K) - beta*R*W == 0;
eq2 = delta*R*W - gamma*W == 0;
% Solve the equations
equilibriumPoints = solve([eq1, eq2], [R, W]);
% Extract the equilibrium point values
Req = double(equilibriumPoints.R);
Weq = double(equilibriumPoints.W);
% Display the equilibrium points
equilibriumPointsValues = [Req, Weq]
% Solve the differential equations using ode45
modified_lv = @(t,Y)[alpha*Y(1)*(1-Y(1)/K)-beta*Y(1)*Y(2);
delta*Y(1)*Y(2)-gamma*Y(2)];
% Initial conditions
R0 = 40;
W0 = 9;
Y0 = [R0, W0];
tspan = [0, 100];
% Solve the differential equations
[t, Y] = ode45(modified_lv, tspan, Y0);
% Extract the populations
R = Y(:, 1);
W = Y(:, 2);
% Plot the results
figure;
subplot(2,1,1);
plot(t, R, 'r', 'LineWidth', 1.5);
hold on;
plot(t, W, 'b', 'LineWidth', 1.5);
xlabel('Time (months)');
ylabel('Population');
legend('R', 'W');
grid on;
subplot(2,1,2);
plot(R, W, 'k', 'LineWidth', 1.5);
xlabel('R');
ylabel('W');
grid on;
hold on;
plot(Req, Weq, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');
legend('Phase Trajectory', 'Equilibrium Point');
Hi everyone,
I've recently joined a forest protection team in Greece, where we use drones for various tasks. This has sparked my interest in drone programming, and I'd like to learn more about it. Can anyone recommend any beginner-friendly courses or programs that teach drone programming?
I'm particularly interested in courses that focus on practical applications and might align with the work we do in forest protection. Any suggestions or guidance would be greatly appreciated!
Thank you!
Gabriel's horn is a shape with the paradoxical property that it has infinite surface area, but a finite volume.
Gabriel’s horn is formed by taking the graph of 
with the domain 
and rotating it in three dimensions about the 
axis.
with the domain and rotating it in three dimensions about the axis.There is a standard formula for calculating the volume of this shape, for a general function 
.Wwe will just state that the volume of the 
solid between a and b is:
.Wwe will just state that the volume of the solid between a and b is:
The surface area of the solid is given by:
One other thing we need to consider is that we are trying to find the value of these integrals between 1 and ∞. An integral with a limit of infinity is called an improper integral and we can't evaluate it simply by plugging the value infinity into the normal equation for a definite integral. Instead, we must first calculate the definite integral up to some finite limit b and then calculate the limit of the result as b tends to ∞:
Volume
We can calculate the horn's volume using the volume integral above, so
The total volume of this infinitely long trumpet isπ.
Surface Area
To determine the surface area, we first need the function’s derivative:
Now plug it into the surface area formula and we have:
This is an improper integral and it's hard to evaluate, but since in our interval 
So, we have :
Now,we evaluate this last integral
So the surface are is infinite.
% Define the function for Gabriel's Horn
gabriels_horn = @(x) 1 ./ x;
% Create a range of x values
x = linspace(1, 40, 4000); % Increase the number of points for better accuracy
y = gabriels_horn(x);
% Create the meshgrid
theta = linspace(0, 2 * pi, 6000); % Increase theta points for a smoother surface
[X, T] = meshgrid(x, theta);
Y = gabriels_horn(X) .* cos(T);
Z = gabriels_horn(X) .* sin(T);
% Plot the surface of Gabriel's Horn
figure('Position', [200, 100, 1200, 900]);
surf(X, Y, Z, 'EdgeColor', 'none', 'FaceAlpha', 0.9);
hold on;
% Plot the central axis
plot3(x, zeros(size(x)), zeros(size(x)), 'r', 'LineWidth', 2);
% Set labels
xlabel('x');
ylabel('y');
zlabel('z');
% Adjust colormap and axis properties
colormap('gray');
shading interp; % Smooth shading
% Adjust the view
view(3);
axis tight;
grid on;
% Add formulas as text annotations
dim1 = [0.4 0.7 0.3 0.2];
annotation('textbox',dim1,'String',{'$$V = \pi \int_{1}^{a} \left( \frac{1}{x} \right)^2 dx = \pi \left( 1 - \frac{1}{a} \right)$$', ...
'', ... % Add an empty line for larger gap
'$$\lim_{a \to \infty} V = \lim_{a \to \infty} \pi \left( 1 - \frac{1}{a} \right) = \pi$$'}, ...
'Interpreter','latex','FontSize',12, 'EdgeColor','none', 'FitBoxToText', 'on');
dim2 = [0.4 0.5 0.3 0.2];
annotation('textbox',dim2,'String',{'$$A = 2\pi \int_{1}^{a} \frac{1}{x} \sqrt{1 + \left( -\frac{1}{x^2} \right)^2} dx > 2\pi \int_{1}^{a} \frac{dx}{x} = 2\pi \ln(a)$$', ...
'', ... % Add an empty line for larger gap
'$$\lim_{a \to \infty} A \geq \lim_{a \to \infty} 2\pi \ln(a) = \infty$$'}, ...
'Interpreter','latex','FontSize',12, 'EdgeColor','none', 'FitBoxToText', 'on');
% Add Gabriel's Horn label
dim3 = [0.3 0.9 0.3 0.1];
annotation('textbox',dim3,'String','Gabriel''s Horn', ...
'Interpreter','latex','FontSize',14, 'EdgeColor','none', 'HorizontalAlignment', 'center');
hold off
daspect([3.5 1 1]) % daspect([x y z])
view(-27, 15)
lightangle(-50,0)
lighting('gouraud')
The properties of this figure were first studied by Italian physicist and mathematician Evangelista Torricelli in the 17th century.
Acknowledgment
I would like to express my sincere gratitude to all those who have supported and inspired me throughout this project.
First and foremost, I would like to thank the mathematician and my esteemed colleague, Stavros Tsalapatis, for inspiring me with the fascinating subject of Gabriel's Horn.
I am also deeply thankful to Mr. @Star Strider for his invaluable assistance in completing the final code.
References:
We are modeling the introduction of a novel pathogen into a completely susceptible population. In the cells below, I have provided you with the Matlab code for a simple stochastic SIR model, implemented using the "GillespieSSA" function
Simulating the stochastic model 100 times for 
Since γ is 0.4 per day, 
per day
per day% Define the parameters
beta = 0.36;
gamma = 0.4;
n_sims = 100;
tf = 100; % Time frame changed to 100
% Calculate R0
R0 = beta / gamma
% Initial state values
initial_state_values = [1000000; 1; 0; 0]; % S, I, R, cum_inc
% Define the propensities and state change matrix
a = @(state) [beta * state(1) * state(2) / 1000000, gamma * state(2)];
nu = [-1, 0; 1, -1; 0, 1; 0, 0];
% Define the Gillespie algorithm function
function [t_values, state_values] = gillespie_ssa(initial_state, a, nu, tf)
t = 0;
state = initial_state(:); % Ensure state is a column vector
t_values = t;
state_values = state';
while t < tf
rates = a(state);
rate_sum = sum(rates);
if rate_sum == 0
break;
end
tau = -log(rand) / rate_sum;
t = t + tau;
r = rand * rate_sum;
cum_sum_rates = cumsum(rates);
reaction_index = find(cum_sum_rates >= r, 1);
state = state + nu(:, reaction_index);
% Update cumulative incidence if infection occurred
if reaction_index == 1
state(4) = state(4) + 1; % Increment cumulative incidence
end
t_values = [t_values; t];
state_values = [state_values; state'];
end
end
% Function to simulate the stochastic model multiple times and plot results
function simulate_stoch_model(beta, gamma, n_sims, tf, initial_state_values, R0, plot_type)
% Define the propensities and state change matrix
a = @(state) [beta * state(1) * state(2) / 1000000, gamma * state(2)];
nu = [-1, 0; 1, -1; 0, 1; 0, 0];
% Set random seed for reproducibility
rng(11);
% Initialize plot
figure;
hold on;
for i = 1:n_sims
[t, output] = gillespie_ssa(initial_state_values, a, nu, tf);
% Check if the simulation had only one step and re-run if necessary
while length(t) == 1
[t, output] = gillespie_ssa(initial_state_values, a, nu, tf);
end
if strcmp(plot_type, 'cumulative_incidence')
plot(t, output(:, 4), 'LineWidth', 2, 'Color', rand(1, 3));
elseif strcmp(plot_type, 'prevalence')
plot(t, output(:, 2), 'LineWidth', 2, 'Color', rand(1, 3));
end
end
xlabel('Time (days)');
if strcmp(plot_type, 'cumulative_incidence')
ylabel('Cumulative Incidence');
ylim([0 inf]);
elseif strcmp(plot_type, 'prevalence')
ylabel('Prevalence of Infection');
ylim([0 50]);
end
title(['Stochastic model output for R0 = ', num2str(R0)]);
subtitle([num2str(n_sims), ' simulations']);
xlim([0 tf]);
grid on;
hold off;
end
% Simulate the model 100 times and plot cumulative incidence
simulate_stoch_model(beta, gamma, n_sims, tf, initial_state_values, R0, 'cumulative_incidence');
% Simulate the model 100 times and plot prevalence
simulate_stoch_model(beta, gamma, n_sims, tf, initial_state_values, R0, 'prevalence');
While searching the internet for some books on ordinary differential equations, I came across a link that I believe is very useful for all math students and not only. If you are interested in ODEs, it's worth taking the time to study it.
A First Look at Ordinary Differential Equations by Timothy S. Judson is an excellent resource for anyone looking to understand ODEs better. Here's a brief overview of the main topics covered:
- Introduction to ODEs: Basic concepts, definitions, and initial differential equations.
- Methods of Solution:
- Separable equations
- First-order linear equations
- Exact equations
- Transcendental functions
- Applications of ODEs: Practical examples and applications in various scientific fields.
- Systems of ODEs: Analysis and solutions of systems of differential equations.
- Series and Numerical Methods: Use of series and numerical methods for solving ODEs.
This book provides a clear and comprehensive introduction to ODEs, making it suitable for students and new researchers in mathematics. If you're interested, you can explore the book in more detail here: A First Look at Ordinary Differential Equations.
The study of the dynamics of the discrete Klein - Gordon equation (DKG) with friction is given by the equation :
In the above equation, W describes the potential function:
to which every coupled unit 
adheres. In Eq. (1), the variable $
$ is the unknown displacement of the oscillator occupying the n-th position of the lattice, and 
is the discretization parameter. We denote by h the distance between the oscillators of the lattice. The chain (DKG) contains linear damping with a damping coefficient 
, while
is the coefficient of the nonlinear cubic term.
adheres. In Eq. (1), the variable $$ is the unknown displacement of the oscillator occupying the n-th position of the lattice, and is the discretization parameter. We denote by h the distance between the oscillators of the lattice. The chain (DKG) contains linear damping with a damping coefficient , whileis the coefficient of the nonlinear cubic term.For the DKG chain (1), we will consider the problem of initial-boundary values, with initial conditions
and Dirichlet boundary conditions at the boundary points 
and 
, that is,
and , that is,
Therefore, when necessary, we will use the short notation 
for the one-dimensional discrete Laplacian
for the one-dimensional discrete Laplacian
Now we want to investigate numerically the dynamics of the system (1)-(2)-(3). Our first aim is to conduct a numerical study of the property of Dynamic Stability of the system, which directly depends on the existence and linear stability of the branches of equilibrium points.
For the discussion of numerical results, it is also important to emphasize the role of the parameter 
. By changing the time variable 
, we rewrite Eq. (1) in the form
. By changing the time variable , we rewrite Eq. (1) in the form
. We consider spatially extended initial conditions of the form:
where 
is the distance of the grid and 
is the amplitude of the initial condition We also assume zero initial velocity:
the following graphs for 
and 
and % Parameters
L = 200; % Length of the system
K = 99; % Number of spatial points
j = 2; % Mode number
omega_d = 1; % Characteristic frequency
beta = 1; % Nonlinearity parameter
delta = 0.05; % Damping coefficient
% Spatial grid
h = L / (K + 1);
n = linspace(-L/2, L/2, K+2); % Spatial points
N = length(n);
omegaDScaled = h * omega_d;
deltaScaled = h * delta;
% Time parameters
dt = 1; % Time step
tmax = 3000; % Maximum time
tspan = 0:dt:tmax; % Time vector
% Values of amplitude 'a' to iterate over
a_values = [2, 1.95, 1.9, 1.85, 1.82]; % Modify this array as needed
% Differential equation solver function
function dYdt = odefun(~, Y, N, h, omegaDScaled, deltaScaled, beta)
U = Y(1:N);
Udot = Y(N+1:end);
Uddot = zeros(size(U));
% Laplacian (discrete second derivative)
for k = 2:N-1
Uddot(k) = (U(k+1) - 2 * U(k) + U(k-1)) ;
end
% System of equations
dUdt = Udot;
dUdotdt = Uddot - deltaScaled * Udot + omegaDScaled^2 * (U - beta * U.^3);
% Pack derivatives
dYdt = [dUdt; dUdotdt];
end
% Create a figure for subplots
figure;
% Initial plot
a_init = 2; % Example initial amplitude for the initial condition plot
U0_init = a_init * sin((j * pi * h * n) / L); % Initial displacement
U0_init(1) = 0; % Boundary condition at n = 0
U0_init(end) = 0; % Boundary condition at n = K+1
subplot(3, 2, 1);
plot(n, U0_init, 'r.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
xlabel('$x_n$', 'Interpreter', 'latex');
ylabel('$U_n$', 'Interpreter', 'latex');
title('$t=0$', 'Interpreter', 'latex');
set(gca, 'FontSize', 12, 'FontName', 'Times');
xlim([-L/2 L/2]);
ylim([-3 3]);
grid on;
% Loop through each value of 'a' and generate the plot
for i = 1:length(a_values)
a = a_values(i);
% Initial conditions
U0 = a * sin((j * pi * h * n) / L); % Initial displacement
U0(1) = 0; % Boundary condition at n = 0
U0(end) = 0; % Boundary condition at n = K+1
Udot0 = zeros(size(U0)); % Initial velocity
% Pack initial conditions
Y0 = [U0, Udot0];
% Solve ODE
opts = odeset('RelTol', 1e-5, 'AbsTol', 1e-6);
[t, Y] = ode45(@(t, Y) odefun(t, Y, N, h, omegaDScaled, deltaScaled, beta), tspan, Y0, opts);
% Extract solutions
U = Y(:, 1:N);
Udot = Y(:, N+1:end);
% Plot final displacement profile
subplot(3, 2, i+1);
plot(n, U(end,:), 'b.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
xlabel('$x_n$', 'Interpreter', 'latex');
ylabel('$U_n$', 'Interpreter', 'latex');
title(['$t=3000$, $a=', num2str(a), '$'], 'Interpreter', 'latex');
set(gca, 'FontSize', 12, 'FontName', 'Times');
xlim([-L/2 L/2]);
ylim([-2 2]);
grid on;
end
% Adjust layout
set(gcf, 'Position', [100, 100, 1200, 900]); % Adjust figure size as needed
Dynamics for the initial condition , 
, for 
, for different amplitude values. By reducing the amplitude values, we observe the convergence to equilibrium points of different branches from 
and the appearance of values 
for which the solution converges to a non-linear equilibrium point 
Parameters: 
, for , for different amplitude values. By reducing the amplitude values, we observe the convergence to equilibrium points of different branches from and the appearance of values for which the solution converges to a non-linear equilibrium point Parameters: 
Detection of a stability threshold 
: For 
, the initial condition , 
, converges to a non-linear equilibrium point
.
: For , the initial condition , , converges to a non-linear equilibrium point.Characteristics for 
, with corresponding norm 
where the dynamics appear in the first image of the third row, we observe convergence to a non-linear equilibrium point of branch 
This has the same norm and the same energy as the previous case but the final state has a completely different profile. This result suggests secondary bifurcations have occurred in branch 
, with corresponding norm where the dynamics appear in the first image of the third row, we observe convergence to a non-linear equilibrium point of branch This has the same norm and the same energy as the previous case but the final state has a completely different profile. This result suggests secondary bifurcations have occurred in branch By further reducing the amplitude, distinct values of 
are discerned: 1.9, 1.85, 1.81 for which the initial condition 
with norms 
respectively, converges to a non-linear equilibrium point of branch 
This equilibrium point has norm 
and energy 
. The behavior of this equilibrium is illustrated in the third row and in the first image of the third row of Figure 1, and also in the first image of the third row of Figure 2. For all the values between the aforementioned a, the initial condition 
converges to geometrically different non-linear states of branch 
as shown in the second image of the first row and the first image of the second row of Figure 2, for amplitudes 
and 
respectively.
are discerned: 1.9, 1.85, 1.81 for which the initial condition with norms respectively, converges to a non-linear equilibrium point of branch This equilibrium point has norm and energy . The behavior of this equilibrium is illustrated in the third row and in the first image of the third row of Figure 1, and also in the first image of the third row of Figure 2. For all the values between the aforementioned a, the initial condition converges to geometrically different non-linear states of branch as shown in the second image of the first row and the first image of the second row of Figure 2, for amplitudes and respectively. Refference:
The study of the dynamics of the discrete Klein - Gordon equation (DKG) with friction is given by the equation :
above equation, W describes the potential function :
The objective of this simulation is to model the dynamics of a segment of DNA under thermal fluctuations with fixed boundaries using a modified discrete Klein-Gordon equation. The model incorporates elasticity, nonlinearity, and damping to provide insights into the mechanical behavior of DNA under various conditions.
% Parameters
numBases = 200; % Number of base pairs, representing a segment of DNA
kappa = 0.1; % Elasticity constant
omegaD = 0.2; % Frequency term
beta = 0.05; % Nonlinearity coefficient
delta = 0.01; % Damping coefficient
- Position: Random initial perturbations between 0.01 and 0.02 to simulate the thermal fluctuations at the start.
- Velocity: All bases start from rest, assuming no initial movement except for the thermal perturbations.
% Random initial perturbations to simulate thermal fluctuations
initialPositions = 0.01 + (0.02-0.01).*rand(numBases,1);
initialVelocities = zeros(numBases,1); % Assuming initial rest state
The simulation uses fixed ends to model the DNA segment being anchored at both ends, which is typical in experimental setups for studying DNA mechanics. The equations of motion for each base are derived from a modified discrete Klein-Gordon equation with the inclusion of damping:
% Define the differential equations
dt = 0.05; % Time step
tmax = 50; % Maximum time
tspan = 0:dt:tmax; % Time vector
x = zeros(numBases, length(tspan)); % Displacement matrix
x(:,1) = initialPositions; % Initial positions
% Velocity-Verlet algorithm for numerical integration
for i = 2:length(tspan)
% Compute acceleration for internal bases
acceleration = zeros(numBases,1);
for n = 2:numBases-1
acceleration(n) = kappa * (x(n+1, i-1) - 2 * x(n, i-1) + x(n-1, i-1)) ...
- delta * initialVelocities(n) - omegaD^2 * (x(n, i-1) - beta * x(n, i-1)^3);
end
% positions for internal bases
x(2:numBases-1, i) = x(2:numBases-1, i-1) + dt * initialVelocities(2:numBases-1) ...
+ 0.5 * dt^2 * acceleration(2:numBases-1);
% velocities using new accelerations
newAcceleration = zeros(numBases,1);
for n = 2:numBases-1
newAcceleration(n) = kappa * (x(n+1, i) - 2 * x(n, i) + x(n-1, i)) ...
- delta * initialVelocities(n) - omegaD^2 * (x(n, i) - beta * x(n, i)^3);
end
initialVelocities(2:numBases-1) = initialVelocities(2:numBases-1) + 0.5 * dt * (acceleration(2:numBases-1) + newAcceleration(2:numBases-1));
end
% Visualization of displacement over time for each base pair
figure;
hold on;
for n = 2:numBases-1
plot(tspan, x(n, :));
end
xlabel('Time');
ylabel('Displacement');
legend(arrayfun(@(n) ['Base ' num2str(n)], 2:numBases-1, 'UniformOutput', false));
title('Displacement of DNA Bases Over Time');
hold off;
The results are visualized using a plot that shows the displacements of each base over time . Key observations from the simulation include :
- Wave Propagation: The initial perturbations lead to wave-like dynamics along the segment, with visible propagation and reflection at the boundaries.
- Damping Effects: The inclusion of damping leads to a gradual reduction in the amplitude of the oscillations, indicating energy dissipation over time.
- Nonlinear Behavior: The nonlinear term influences the response, potentially stabilizing the system against large displacements or leading to complex dynamic patterns.
% 3D plot for displacement
figure;
[X, T] = meshgrid(1:numBases, tspan);
surf(X', T', x);
xlabel('Base Pair');
ylabel('Time');
zlabel('Displacement');
title('3D View of DNA Base Displacements');
colormap('jet');
shading interp;
colorbar; % Adds a color bar to indicate displacement magnitude
% Snapshot visualization at a specific time
snapshotTime = 40; % Desired time for the snapshot
[~, snapshotIndex] = min(abs(tspan - snapshotTime)); % Find closest index
snapshotSolution = x(:, snapshotIndex); % Extract displacement at the snapshot time
% Plotting the snapshot
figure;
stem(1:numBases, snapshotSolution, 'filled'); % Discrete plot using stem
title(sprintf('DNA Model Displacement at t = %d seconds', snapshotTime));
xlabel('Base Pair Index');
ylabel('Displacement');
% Time vector for detailed sampling
tDetailed = 0:0.5:50; % Detailed time steps
% Initialize an empty array to hold the data
data = [];
% Generate the data for 3D plotting
for i = 1:numBases
% Interpolate to get detailed solution data for each base pair
detailedSolution = interp1(tspan, x(i, :), tDetailed);
% Concatenate the current base pair's data to the main data array
data = [data; repmat(i, length(tDetailed), 1), tDetailed', detailedSolution'];
end
% 3D Plot
figure;
scatter3(data(:,1), data(:,2), data(:,3), 10, data(:,3), 'filled');
xlabel('Base Pair');
ylabel('Time');
zlabel('Displacement');
title('3D Plot of DNA Base Pair Displacements Over Time');
colorbar; % Adds a color bar to indicate displacement magnitude
Hello MathWorks Community,
I am excited to announce that I am currently working on a book project centered around Matrix Algebra, specifically designed for MATLAB users. This book aims to cater to undergraduate students in engineering, where Matrix Algebra serves as a foundational element.
Matrix Algebra is not only pivotal in understanding complex engineering concepts but also in applying these principles effectively in various technological solutions. MATLAB, renowned for its powerful computational capabilities, is an excellent tool to explore and implement these concepts, making it a perfect companion for this book.
As I embark on this journey to create a resource that bridges theoretical matrix algebra with practical MATLAB applications, I am looking for one or two knowledgeable individuals who have a firm grasp of both subjects. If you have experience in teaching or applying matrix algebra in engineering contexts and are familiar with MATLAB, your contribution could be invaluable.
Collaborators will help in shaping the content to ensure it is educational, engaging, and technically robust, making complex concepts accessible and applicable for students.
If you are interested in contributing to this project or know someone who might be, please reach out to discuss how we can work together to make this book a valuable resource for engineering students.
Thank you and looking forward to your participation!
What's your way?
I'm excited to share some valuable resources that I've found to be incredibly helpful for anyone looking to enhance their MATLAB skills. Whether you're just starting out, studying as a student, or are a seasoned professional, these guides and books offer a wealth of information to aid in your learning journey.
These materials are freely available and can be a great addition to your learning resources. They cover a wide range of topics and are designed to help users at all levels to improve their proficiency in MATLAB.
Happy learning and I hope you find these resources as useful as I have!
Let S be the closed surface composed of the hemisphere 
and the base 
Let 
be the electric field defined by 
. Find the electric flux through S. (Hint: Divide S into two parts and calculate 
).
and the base Let be the electric field defined by . Find the electric flux through S. (Hint: Divide S into two parts and calculate ).% Define the limits of integration for the hemisphere S1
theta_lim = [-pi/2, pi/2];
phi_lim = [0, pi/2];
% Perform the double integration over the spherical surface of the hemisphere S1
% Define the electric flux function for the hemisphere S1
flux_function_S1 = @(theta, phi) 2 * sin(phi);
electric_flux_S1 = integral2(flux_function_S1, theta_lim(1), theta_lim(2), phi_lim(1), phi_lim(2));
% For the base of the hemisphere S2, the electric flux is 0 since the electric
% field has no z-component at the base
electric_flux_S2 = 0;
% Calculate the total electric flux through the closed surface S
total_electric_flux = electric_flux_S1 + electric_flux_S2;
% Display the flux calculations
disp(['Electric flux through the hemisphere S1: ', num2str(electric_flux_S1)]);
disp(['Electric flux through the base of the hemisphere S2: ', num2str(electric_flux_S2)]);
disp(['Total electric flux through the closed surface S: ', num2str(total_electric_flux)]);
% Parameters for the plot
radius = 1; % Radius of the hemisphere
% Create a meshgrid for theta and phi for the plot
[theta, phi] = meshgrid(linspace(theta_lim(1), theta_lim(2), 20), linspace(phi_lim(1), phi_lim(2), 20));
% Calculate Cartesian coordinates for the points on the hemisphere
x = radius * sin(phi) .* cos(theta);
y = radius * sin(phi) .* sin(theta);
z = radius * cos(phi);
% Define the electric field components
Ex = 2 * x;
Ey = 2 * y;
Ez = 2 * z;
% Plot the hemisphere
figure;
surf(x, y, z, 'FaceAlpha', 0.5, 'EdgeColor', 'none');
hold on;
% Plot the electric field vectors
quiver3(x, y, z, Ex, Ey, Ez, 'r');
% Plot the base of the hemisphere
[x_base, y_base] = meshgrid(linspace(-radius, radius, 20), linspace(-radius, radius, 20));
z_base = zeros(size(x_base));
surf(x_base, y_base, z_base, 'FaceColor', 'cyan', 'FaceAlpha', 0.3);
% Additional plot settings
colormap('cool');
axis equal;
grid on;
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Hemisphere and Electric Field');
The line integral 
, where C is the boundary of the square 
oriented counterclockwise, can be evaluated in two ways:
, where C is the boundary of the square oriented counterclockwise, can be evaluated in two ways: Using the definition of the line integral:
% Initialize the integral sum
integral_sum = 0;
% Segment C1: x = -1, y goes from -1 to 1
y = linspace(-1, 1);
x = -1 * ones(size(y));
dy = diff(y);
integral_sum = integral_sum + sum(-x(1:end-1) .* dy);
% Segment C2: y = 1, x goes from -1 to 1
x = linspace(-1, 1);
y = ones(size(x));
dx = diff(x);
integral_sum = integral_sum + sum(y(1:end-1).^2 .* dx);
% Segment C3: x = 1, y goes from 1 to -1
y = linspace(1, -1);
x = ones(size(y));
dy = diff(y);
integral_sum = integral_sum + sum(-x(1:end-1) .* dy);
% Segment C4: y = -1, x goes from 1 to -1
x = linspace(1, -1);
y = -1 * ones(size(x));
dx = diff(x);
integral_sum = integral_sum + sum(y(1:end-1).^2 .* dx);
disp(['Direct Method Integral: ', num2str(integral_sum)]);
Plotting the square path
% Define the square's vertices
vertices = [-1 -1; -1 1; 1 1; 1 -1; -1 -1];
% Plot the square
figure;
plot(vertices(:,1), vertices(:,2), '-o');
title('Square Path for Line Integral');
xlabel('x');
ylabel('y');
grid on;
axis equal;
% Add arrows to indicate the path direction (counterclockwise)
hold on;
for i = 1:size(vertices,1)-1
% Calculate direction
dx = vertices(i+1,1) - vertices(i,1);
dy = vertices(i+1,2) - vertices(i,2);
% Reduce the length of the arrow for better visibility
scale = 0.2;
dx = scale * dx;
dy = scale * dy;
% Calculate the start point of the arrow
startx = vertices(i,1) + (1 - scale) * dx;
starty = vertices(i,2) + (1 - scale) * dy;
% Plot the arrow
quiver(startx, starty, dx, dy, 'MaxHeadSize', 0.5, 'Color', 'r', 'AutoScale', 'off');
end
hold off;
Apply Green's Theorem for the line integral
% Define the partial derivatives of P and Q
f = @(x, y) -1 - 2*y; % derivative of -x with respect to x is -1, and derivative of y^2 with respect to y is 2y
% Compute the double integral over the square [-1,1]x[-1,1]
integral_value = integral2(f, -1, 1, 1, -1);
disp(['Green''s Theorem Integral: ', num2str(integral_value)]);
Plotting the vector field related to Green’s theorem
% Define the grid for the vector field
[x, y] = meshgrid(linspace(-2, 2, 20), linspace(-2 ,2, 20));
% Define the vector field components
P = y.^2; % y^2 component
Q = -x; % -x component
% Plot the vector field
figure;
quiver(x, y, P, Q, 'b');
hold on; % Hold on to plot the square on the same figure
% Define the square's vertices
vertices = [-1 -1; -1 1; 1 1; 1 -1; -1 -1];
% Plot the square path
plot(vertices(:,1), vertices(:,2), '-o', 'Color', 'k'); % 'k' for black color
title('Vector Field (P = y^2, Q = -x) with Square Path');
xlabel('x');
ylabel('y');
axis equal;
% Add arrows to indicate the path direction (counterclockwise)
for i = 1:size(vertices,1)-1
% Calculate direction
dx = vertices(i+1,1) - vertices(i,1);
dy = vertices(i+1,2) - vertices(i,2);
% Reduce the length of the arrow for better visibility
scale = 0.2;
dx = scale * dx;
dy = scale * dy;
% Calculate the start point of the arrow
startx = vertices(i,1) + (1 - scale) * dx;
starty = vertices(i,2) + (1 - scale) * dy;
% Plot the arrow
quiver(startx, starty, dx, dy, 'MaxHeadSize', 0.5, 'Color', 'r', 'AutoScale', 'off');
end
hold off;
To solve a surface integral for example the
over the sphere 
easily in MATLAB, you can leverage the symbolic toolbox for a direct and clear solution. Here is a tip to simplify the process:
over the sphere easily in MATLAB, you can leverage the symbolic toolbox for a direct and clear solution. Here is a tip to simplify the process:- Use Symbolic Variables and Functions: Define your variables symbolically, including the parameters of your spherical coordinates θ and ϕ and the radius r . This allows MATLAB to handle the expressions symbolically, making it easier to manipulate and integrate them.
- Express in Spherical Coordinates Directly: Since you already know the sphere's equation and the relationship in spherical coordinates, define x, y, and z in terms of r , θ and ϕ directly.
- Perform Symbolic Integration: Use MATLAB's `int` function to integrate symbolically. Since the sphere and the function
are symmetric, you can exploit these symmetries to simplify the calculation.
Here’s how you can apply this tip in MATLAB code:
% Include the symbolic math toolbox
syms theta phi
% Define the limits for theta and phi
theta_limits = [0, pi];
phi_limits = [0, 2*pi];
% Define the integrand function symbolically
integrand = 16 * sin(theta)^3 * cos(phi)^2;
% Perform the symbolic integral for the surface integral
surface_integral = int(int(integrand, theta, theta_limits(1), theta_limits(2)), phi, phi_limits(1), phi_limits(2));
% Display the result of the surface integral symbolically
disp(['The surface integral of x^2 over the sphere is ', char(surface_integral)]);
% Number of points for plotting
num_points = 100;
% Define theta and phi for the sphere's surface
[theta_mesh, phi_mesh] = meshgrid(linspace(double(theta_limits(1)), double(theta_limits(2)), num_points), ...
linspace(double(phi_limits(1)), double(phi_limits(2)), num_points));
% Spherical to Cartesian conversion for plotting
r = 2; % radius of the sphere
x = r * sin(theta_mesh) .* cos(phi_mesh);
y = r * sin(theta_mesh) .* sin(phi_mesh);
z = r * cos(theta_mesh);
% Plot the sphere
figure;
surf(x, y, z, 'FaceColor', 'interp', 'EdgeColor', 'none');
colormap('jet'); % Color scheme
shading interp; % Smooth shading
camlight headlight; % Add headlight-type lighting
lighting gouraud; % Use Gouraud shading for smooth color transitions
title('Sphere: x^2 + y^2 + z^2 = 4');
xlabel('x-axis');
ylabel('y-axis');
zlabel('z-axis');
colorbar; % Add color bar to indicate height values
axis square; % Maintain aspect ratio to be square
view([-30, 20]); % Set a nice viewing angle
Before we begin, you will need to make sure you have 'sir_age_model.m' installed. Once you've downloaded this folder into your working directory, which can be located at your current folder. If you can see this file in your current folder, then it's safe to use it. If you choose to use MATLAB online or MATLAB Mobile, you may upload this to your MATLAB Drive.
This is the code for the SIR model stratified into 2 age groups (children and adults). For a detailed explanation of how to derive the force of infection by age group.
% Main script to run the SIR model simulation
% Initial state values
initial_state_values = [200000; 1; 0; 800000; 0; 0]; % [S1; I1; R1; S2; I2; R2]
% Parameters
parameters = [0.05; 7; 6; 1; 10; 1/5]; % [b; c_11; c_12; c_21; c_22; gamma]
% Time span for the simulation (3 months, with daily steps)
tspan = [0 90];
% Solve the ODE
[t, y] = ode45(@(t, y) sir_age_model(t, y, parameters), tspan, initial_state_values);
% Plotting the results
plot(t, y);
xlabel('Time (days)');
ylabel('Number of people');
legend('S1', 'I1', 'R1', 'S2', 'I2', 'R2');
title('SIR Model with Age Structure');
What was the cumulative incidence of infection during this epidemic? What proportion of those infections occurred in children?
In the SIR model, the cumulative incidence of infection is simply the decline in susceptibility.
% Assuming 'y' contains the simulation results from the ode45 function
% and 't' contains the time points
% Total cumulative incidence
total_cumulative_incidence = (y(1,1) - y(end,1)) + (y(1,4) - y(end,4));
fprintf('Total cumulative incidence: %f\n', total_cumulative_incidence);
% Cumulative incidence in children
cumulative_incidence_children = (y(1,1) - y(end,1));
% Proportion of infections in children
proportion_infections_children = cumulative_incidence_children / total_cumulative_incidence;
fprintf('Proportion of infections in children: %f\n', proportion_infections_children);
927,447 people became infected during this epidemic, 20.5% of which were children.
Which age group was most affected by the epidemic?
To answer this, we can calculate the proportion of children and adults that became infected.
% Assuming 'y' contains the simulation results from the ode45 function
% and 't' contains the time points
% Proportion of children that became infected
initial_children = 200000; % initial number of susceptible children
final_susceptible_children = y(end,1); % final number of susceptible children
proportion_infected_children = (initial_children - final_susceptible_children) / initial_children;
fprintf('Proportion of children that became infected: %f\n', proportion_infected_children);
% Proportion of adults that became infected
initial_adults = 800000; % initial number of susceptible adults
final_susceptible_adults = y(end,4); % final number of susceptible adults
proportion_infected_adults = (initial_adults - final_susceptible_adults) / initial_adults;
fprintf('Proportion of adults that became infected: %f\n', proportion_infected_adults);
Throughout this epidemic, 95% of all children and 92% of all adults were infected. Children were therefore slightly more affected in proportion to their population size, even though the majority of infections occurred in adults.
I would like to propose the creation of MATLAB EduHub, a dedicated channel within the MathWorks community where educators, students, and professionals can share and access a wealth of educational material that utilizes MATLAB. This platform would act as a central repository for articles, teaching notes, and interactive learning modules that integrate MATLAB into the teaching and learning of various scientific fields.
Key Features:
1. Resource Sharing: Users will be able to upload and share their own educational materials, such as articles, tutorials, code snippets, and datasets.
2. Categorization and Search: Materials can be categorized for easy searching by subject area, difficulty level, and MATLAB version..
3. Community Engagement: Features for comments, ratings, and discussions to encourage community interaction.
4. Support for Educators: Special sections for educators to share teaching materials and track engagement.
Benefits:
- Enhanced Educational Experience: The platform will enrich the learning experience through access to quality materials.
- Collaboration and Networking: It will promote collaboration and networking within the MATLAB community.
- Accessibility of Resources: It will make educational materials available to a wider audience.
By establishing MATLAB EduHub, I propose a space where knowledge and experience can be freely shared, enhancing the educational process and the MATLAB community as a whole.
Keep calm and study PDEs
The stationary solutions of the Klein-Gordon equation refer to solutions that are time-independent, meaning they remain constant over time. For the non-linear Klein-Gordon equation you are discussing:
Stationary solutions arise when the time derivative term, 
, is zero, meaning the motion of the system does not change over time. This leads to a static differential equation:
, is zero, meaning the motion of the system does not change over time. This leads to a static differential equation:
This equation describes how particles in the lattice interact with each other and how non-linearity affects the steady state of the system.
The solutions to this equation correspond to the various possible stable equilibrium states of the system, where each represents different static distribution patterns of displacements 
. The specific form of these stationary solutions depends on the system parameters, such as κ , ω, and β , as well as the initial and boundary conditions of the problem.
. The specific form of these stationary solutions depends on the system parameters, such as κ , ω, and β , as well as the initial and boundary conditions of the problem.To find these solutions in a more specific form, one might need to solve the equation using analytical or numerical methods, considering the different cases that could arise in such a non-linear system.
By interpreting the equation in this way, we can relate the dynamics described by the discrete Klein - Gordon equation to the behavior of DNA molecules within a biological system . This analogy allows us to understand the behavior of DNA in terms of concepts from physics and mathematical modeling .
% Parameters
numBases = 100; % Number of spatial points
omegaD = 0.2; % Common parameter for the equation
% Preallocate the array for the function handles
equations = cell(numBases, 1);
% Initial guess for the solution
initialGuess = 0.01 * ones(numBases, 1);
% Parameter sets for kappa and beta
paramSets = [0.1, 0.05; 0.5, 0.05; 0.1, 0.2];
% Prepare figure for subplot
figure;
set(gcf, 'Position', [100, 100, 1200, 400]); % Set figure size
% Newton-Raphson method parameters
maxIterations = 1000;
tolerance = 1e-10;
% Set options for fsolve to use the 'levenberg-marquardt' algorithm
options = optimoptions('fsolve', 'Algorithm', 'levenberg-marquardt', 'MaxIterations', maxIterations, 'FunctionTolerance', tolerance);
for i = 1:size(paramSets, 1)
kappa = paramSets(i, 1);
beta = paramSets(i, 2);
% Define the equations using a function
for n = 2:numBases-1
equations{n} = @(x) -kappa * (x(n+1) - 2 * x(n) + x(n-1)) - omegaD^2 * (x(n) - beta * x(n)^3);
end
% Boundary conditions with specified fixed values
someFixedValue1 = 10; % Replace with actual value if needed
someFixedValue2 = 10; % Replace with actual value if needed
equations{1} = @(x) x(1) - someFixedValue1;
equations{numBases} = @(x) x(numBases) - someFixedValue2;
% Combine all equations into a single function
F = @(x) cell2mat(cellfun(@(f) f(x), equations, 'UniformOutput', false));
% Solve the system of equations using fsolve with the specified options
x_solution = fsolve(F, initialGuess, options);
norm(F(x_solution))
% Plot the solution in a subplot
subplot(1, 3, i);
plot(x_solution, 'o-', 'LineWidth', 2);
grid on;
xlabel('n', 'FontSize', 12);
ylabel('x[n]', 'FontSize', 12);
title(sprintf('\\kappa = %.2f, \\beta = %.2f', kappa, beta), 'FontSize', 14);
end
% Improve overall aesthetics
sgtitle('Stationary States for Different \kappa and \beta Values', 'FontSize', 16); % Super title for the figure
In the second plot, the elasticity constant κis increased to 0.5, representing a system with greater stiffness . This parameter influences how resistant the system is to deformation, implying that a higher κ makes the system more resilient to changes . By increasing κ, we are essentially tightening the interactions between adjacent units in the model, which could represent, for instance, stronger bonding forces in a physical or biological system .
In the third plot the nonlinearity coefficient β is increased to 0.2 . This adjustment enhances the nonlinear interactions within the system, which can lead to more complex dynamic behaviors, especially in systems exhibiting bifurcations or chaos under certain conditions .
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)