% This file solves the FP equation for an OU process equation


clear
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

L = 4; % length of the domain x in [-L,L]
d = 1; % diffusion coefficient (wlog)
K = 1; % restoring force parameter (wlog)

N=1601;    % number of spatial grid points, an odd number
M = 30; % number of time steps
kstep = 75; %number of intermediate time steps
h =L/((N-1)/2); %grid size
 
Yamp = 0; % this is the amplitude of the forced center y_0(t)
X = h*((1:N)'-(N+1)/2);% grid points

Xm = [X-h/2;X(N)+h/2]; % N+1 midpoints

vel = -K*Xm;
vmx = max(abs(vel));
 
dt = 0.9*h/vmx;  % set time step 
scal = d*dt/h^2;
 % specify initial data

 X0 = 2; % initial position shift
 U0 =  exp(-(X-X0).^2/0.1) ; % a gaussian distribution
 U = U0;
 
 figure(1)
 plot(X,U,'--','linewidth',2 )
 hold on
 
t=0;  % initial time
% preallocate storage
T = zeros(M,1);
inn = zeros(M,1);
mx = zeros(M,1);
mj = zeros(M,1);
 
% uses Crank Nicolson to solve the diffusion equation

% set up the CN matrix

Atm=(1+scal)*ones(N,1);
Atm(1,1) = 1+scal/2;
Atm(N,1) = 1+scal/2;

Am=diag(Atm)-diag(scal/2*ones(N-1,1),1)-diag(scal/2*ones(N-1,1),-1);

 
% The algorithm
inn(1) = sum(U)*h;  % check the total integral
[mx(1),mj(1)] = max(U); % locate the maximum
for n=2:M
for k = 1:kstep
vel = -K*(Xm-Yamp*sin(t/2));
 
%vel = -K*ones(N+1,1);
msk = (vel>0);
% F is the flux
F =  vel.*(msk.*[0;U]+(1-msk).*[U;0]) ; % this uses upwinding
S=0;% source term
Fdiff = (F(2:N+1)-F(1:N))/h;
U_t =  scal/2 * ([U(2:N );U(N )] - 2*U + [U(1 );U(1:N-1 )])-dt*Fdiff+dt*S ;
%U_t =  -dt*Fdiff ;
 
 
U = Am\(U+U_t);
 t=t + dt;
end
T(n) = t;
inn(n) = sum(U)*h;  % check the total integral
[mx(n),mj(n)] = max(U);  % locate the maximum
figure(1)
  plot(X,U,'linewidth',2 )
  %axis([-L/2 L  0 1])
 % axis([0 L -.5 1])
 xlabel('\xi','fontsize',20)
 ylabel('p(\xi,\tau)','fontsize',20)
 hold on
set(gca, 'Xtick',-4:1:4)
end
  
figure(1)
hold off
if(cmykflg==1)
  print('../../figs_c/chapt_10/OU_simulated','-deps','-cmyk')
 end 

figure(2)
plot( T,mx)
ylabel('maximum p(\xi,\tau)','fontsize',20)
figure(3)
plot( T,(X(mj)))
 ylabel('mean(p(\xi,\tau))','fontsize',20)
  