The function delG takes the number of excitations (ne), time step of discretization (delta), and the parameter p (page 354) as input and generates the discretized version of the white noise excitation.
% Code starts from here
function [dw,Ij0,I0j,Ip,I]=delG(ne,p,delta)
% GENERAL CASE
% Ij0=I(j,0) : dz
% I0j=I(0,j)
% Ip=I(j1,j2) : j1,j2=1,2,...ne
% I=I(j1,j2,j3) : j1,j2,j3=1,2,...ne
zii=zeros(ne,1);
si=randn(ne,1);
zi=randn(ne,p);
zi(:,p+1:2*p)=zeros(ne,p);
etajp=randn(ne,p);
etajp(:,p+1:2*p)=zeros(ne,p);
mujp=randn(ne,p);
phyjp=randn(ne,p);
Ap=zeros(ne);
Bp=zeros(ne);
Cp=zeros(ne);
Dp=zeros(ne,ne,ne);
Ip=zeros(ne);
I=zeros(ne,ne,ne);
Jp0=zeros(ne);
Jp=zeros(ne);
J=zeros(ne,ne,ne);
b=zeros(ne,1);
dw=sqrt(delta)*si; % I(j)
r2=0;
r4=0;
for r=1:p
zii=zii+1/r*zi(:,r);
r2=r2+1/r^2;
r4=r4+1/r^4;
end
rhop=(1/12)-(1/(2*pi^2))*r2;
aj0=(-sqrt(2*delta)/pi)*zii-2*sqrt(delta*rhop)*mujp(:,p);
Ij0=delta/2*(sqrt(delta)*si+aj0);
I0j=dw*delta-Ij0;
Ijj=1/2*(dw.^2-delta);
Ijjj=(1/2)*((1/3)*dw.^2-delta).*dw;
alphap=pi^2/180-r4/(2*pi^2);
for j=1:ne
esum=0;
for r=1:p
esum=etajp(j,r)/r^2;
end
b(j)=sqrt(delta/2)*esum+sqrt(delta*alphap)*phyjp(j,p);
end
for j1=1:ne
for j2=1:ne
rsum1=0;
rsum2=0;
rsum3=0;
for r=1:p
rsum1=rsum1+(zi(j1,r)*etajp(j2,r)-zi(j2,r)*etajp(j1,r))/r;
rsum2=(zi(j1,r)*zi(j2,r)+etajp(j1,r)*etajp(j2,r))/r^2;
end
for r=1:p
for l=1:p
if r~=l
rsum3=rsum3+(r/(r^2-l^2))*((1/l)*zi(j1,r)*zi(j2,l)-(l/r)*etajp(j1,r)*etajp(j2,l));
end
end
end
Ap(j1,j2)=1/(2*pi) *rsum1;
Bp(j1,j2)=1/(4*pi^2) *rsum2;
Cp(j1,j2)=-1/(2*pi^2) *rsum3;
Jp(j1,j2)=(1/2)*delta*si(j1)*si(j2)-(1/2)*sqrt(delta)*(aj0(j2)*si(j1)-aj0(j1)*si(j2))+delta*Ap(j1,j2);
Jp0(j1,j2)=(1/factorial(3))*delta^2*si(j1)*si(j2)-delta^1.5/pi*si(j2)*b(j1)+delta^2*Bp(j1,j2)-delta^1.5/4*aj0(j2)*si(j1)+delta^1.5/(2*pi)*si(j1)*b(j2)+delta^2*Cp(j1,j2)+delta^2/2*Ap(j1,j2);
if j1==j2
Ip(j1,j2)=(1/2)*(dw(j1)^2-delta);
else
Ip(j1,j2)=Jp(j1,j2);
end
end
end
for j1=1:ne
for j2=1:ne
for j3=1:ne
if j1==j2
is1=1;
else
is1=0;
end
if j2==j3
is2=1;
else
is2=0;
end
sum1=0;
for r=1:p
for l=1:p
if r<p
t1=zi(j2,l) * (zi(j3,l+r) * etajp(j1,r) - zi(j1,r) * etajp(j1,l+r));
t2=etajp(j2,l) * (etajp(j3,l+r) * etajp(j1,r) + zi(j1,r) * zi(j3,l+r));
sum1=sum1+(1/(l*(l+r))*(t1+t2));
end
end
end
clear t1 t2
term1=sum1;
sum2=0;
for l=2:p
for r=1:(l-1)
if r<p
t1=zi(j2,l)*(zi(j1,r)*etajp(j3,l-r)+zi(j3,l-r)*etajp(j1,r));
t2=etajp(j2,l)*(zi(j1,r)*zi(j3,l-r)-etajp(j1,r)*etajp(j3,l-r));
sum2=sum2+(1/(r*(l-r))*(t1-t2));
end
end
end
clear t1 t2
term2=sum2;
sum3=0;
for l=1:p
for r=(l+1):2*p
if r<p
t1=zi(j2,l)*(zi(j3,r-l)*etajp(j1,r)-zi(j1,r)*etajp(j3,r-l));
t2=etajp(j2,l)*(zi(j1,r)*zi(j3,r-l)+etajp(j1,r)*etajp(j3,r-l));
sum3=sum3+(1/(r*(r-l))*(t1+t2));
end
end
end
term3=sum3;
Dp(j1,j2,j3)=1/(pi^2*2^2.5)*(-term1+term2+term3);
J(j1,j2,j3)=1/sqrt(delta)*si(j1)*Jp0(j2,j3)+1/2*aj0(j1)*Jp(j2,j3)+1/(2*pi)*b(j1)*si(j2)*si(j3)-delta^1.5*si(j2)*Bp(j1,j3)+delta^1.5*si(j3)*(1/2*Ap(j1,j2)-Cp(j2,j1))+delta^1.5*Dp(j1,j2,j3);
I(j1,j2,j3)=J(j1,j2,j3)-0.5*(is1*I0j(j3)+is2*Ij0(j1));
end
end
end
for i=1:ne
I(i,i,i)=Ijjj(i);
end