# The Stochastic Pendulum with red noise

This script is a tutorial script for stochastic integration and visualization using red noise realized using the method of Bartosch.

## The governing ODEs of the damped pendulum are

is taken to vary as where the perturbation is red noise this means no Ito calculus comes in.

## This section is all the book keeping for initialization

Initialize the random number generator and seed it.

rng('shuffle');


Specify the global parameters for the pendulum.

length=1;
g0=9.81;
damp=1e-1;


Specify the time step and stochastic parameters.

numsteps=10; % sets how many time steps between saves.
numouts=500; % sets the number of outputs.
dt=1e-3; % sets the time step.  Stochastic methods are low order so small time steps are required.
numens=1000; % sets the number of ensemble members.


## The parameters for the red noise as realized by Bartosch's algorithm

memt=0.1; % sets the memory of the noise.  Need at least ten time steps in memt.
rho=exp(-dt/memt); rhoc=sqrt(1-rho*rho); % parameters for Bartosch's method.
mysig=0.35; % set the standard deviation of the noise.
prevrand=zeros(numens,1); % start g at g0.


## Initial conditions

%theta=2*pi*rand(numens,1); %initially assume random angles.
theta=(1-1e-4)*pi*ones(numens,1); % start near the unstable state.
%theta=1e-4*pi*ones(numens,1); % start near the stable state.
v=zeros(numens,1); % initially assume zero velocity.
t=0;


## Set up arrays for storage and store initial state

thetas=zeros(numouts+1,numens);
vs=zeros(numouts+1,numens);
ts=zeros(numouts+1,1);
thetas(1,:)=theta;
vs(1,:)=v;
ts(1)=t;


## The main work horse loops. The double loop can be used view to intermediate results

for ii=1:numouts
for jj=1:numsteps
nowrand=randn(numens,1)*mysig; % generate a new bunch of random numbers.
pertnow=(nowrand*rhoc+rho*prevrand); prevrand=pertnow; % use Bartosch's memory to get the right "memory".
g=g0*(1+pertnow); % define g.
% now define the RHS for the DE
rhth=v;
rhv=-(g/length).*sin(theta)-damp*v;
% Step the DE forward. If g had a white noise component this would
% get trickier to account for the sqrt(dt) in Ito's calculus.
theta=theta+dt*rhth;
v=v+dt*rhv;
t=t+dt;
end
thetas(ii+1,:)=theta;
vs(ii+1,:)=v;
ts(ii+1)=t;
end
plot(thetas)


## Now that we have the data we can make a variety of pictures

figure(1)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',14,...
'DefaultTextFontWeight','bold','DefaultAxesFontSize',14,...
'DefaultAxesFontWeight','bold');
subplot(2,1,1)
plot(ts,thetas(:,1))
xlabel('t')
ylabel('theta')
title('one realization')
subplot(2,1,2)
plot(ts,vs(:,1))
xlabel('t')
ylabel('v')

figure(2)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',14,...
'DefaultTextFontWeight','bold','DefaultAxesFontSize',14,...
'DefaultAxesFontWeight','bold');
subplot(2,1,1)
plot(thetas(:,1),vs(:,1))
ylabel('v')
xlabel('theta')
title('one phase space plot')
subplot(2,1,2)
plot(mean(thetas,2),mean(vs,2))
ylabel('v')
xlabel('theta')
title('mean phase space plot')

figure(3)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',14,...
'DefaultTextFontWeight','bold','DefaultAxesFontSize',14,...
'DefaultAxesFontWeight','bold');
subplot(2,1,1)
plot(ts,std(thetas,[],2))
ylabel('standard deviation of theta')
xlabel('t')
subplot(2,1,2)
plot(ts,std(vs,[],2))
ylabel('standard deviation of v')
xlabel('t')

figure(4)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',14,...
'DefaultTextFontWeight','bold','DefaultAxesFontSize',14,...
'DefaultAxesFontWeight','bold');
subplot(2,1,1)
mnth=min(thetas(end,:));
mxth=max(thetas(end,:));
bins=linspace(mnth,mxth,floor(numens/20)+1);
hist(thetas(end,:),bins)
title('histogram of thetas final time')
subplot(2,1,2)
mnv=min(vs(end,:));
mxv=max(vs(end,:));
bins=linspace(mnv,mxv,floor(numens/20)+1);
hist(vs(end,:),bins)
title('histogram of vs final time')