From 57a1f9fb992c090eefd5ace6300fbb662b156958 Mon Sep 17 00:00:00 2001 From: "Nolte, David D" Date: Fri, 29 Jan 2021 11:48:17 -0500 Subject: [PATCH] Add files via upload --- DoublePend.m | 162 +++++++++++++++++++++++++++++++++ Heiles.m | 161 +++++++++++++++++++++++++++++++++ NaK.m | 89 ++++++++++++++++++ PredPrey.m | 47 ++++++++++ avgeigs.m | 50 +++++++++++ continued.m | 20 +++++ convergent.m | 40 +++++++++ coupledrossler.m | 82 +++++++++++++++++ diffusiondriver.m | 158 ++++++++++++++++++++++++++++++++ flowsimple.m | 108 ++++++++++++++++++++++ logistic.m | 40 +++++++++ lorenz.m | 110 +++++++++++++++++++++++ lotkavoltstream.m | 97 ++++++++++++++++++++ quasiHam.m | 145 ++++++++++++++++++++++++++++++ quasiSpec.m | 223 ++++++++++++++++++++++++++++++++++++++++++++++ randsgn.m | 10 +++ repeq.m | 194 ++++++++++++++++++++++++++++++++++++++++ repmut.m | 164 ++++++++++++++++++++++++++++++++++ rksimple.m | 53 +++++++++++ rksimpleVirus2.m | 47 ++++++++++ rossler.m | 93 +++++++++++++++++++ sinecircle1.m | 77 ++++++++++++++++ zerodiag.m | 24 +++++ 23 files changed, 2194 insertions(+) create mode 100644 DoublePend.m create mode 100644 Heiles.m create mode 100644 NaK.m create mode 100644 PredPrey.m create mode 100644 avgeigs.m create mode 100644 continued.m create mode 100644 convergent.m create mode 100644 coupledrossler.m create mode 100644 diffusiondriver.m create mode 100644 flowsimple.m create mode 100644 logistic.m create mode 100644 lorenz.m create mode 100644 lotkavoltstream.m create mode 100644 quasiHam.m create mode 100644 quasiSpec.m create mode 100644 randsgn.m create mode 100644 repeq.m create mode 100644 repmut.m create mode 100644 rksimple.m create mode 100644 rksimpleVirus2.m create mode 100644 rossler.m create mode 100644 sinecircle1.m create mode 100644 zerodiag.m diff --git a/DoublePend.m b/DoublePend.m new file mode 100644 index 0000000..c48a05a --- /dev/null +++ b/DoublePend.m @@ -0,0 +1,162 @@ +function DoublePend +format compact +c(1) = 'r'; +c(2) = 'b'; +c(3) = 'g'; +c(4) = 'k'; +c(5) = 'c'; +c(6) = 'm'; +c(7) = 'y'; + +h = colormap(lines); + +E = 0.8; % 1.0 1.2; % Kinetic Energy + +h2 = figure(2); +h2.Position = [11 -31 560 420]; +dum = set(h2); + +h3 = figure(3); +h3.Position = [10 381 560 420]; +dum = set(h3); + +h5 = figure(5); +h5.Position = [229 1 1008 804]; +dum = set(h5); +clf +hold on + + +repnum = 128; % 32 +mulnum = 64/repnum; +for reploop = 1:repnum + + clear xvar yvar vy vx px py + + vx0 = 2*(0.5 - rand)*sqrt(E); + vy0 = -vx0 + sign(randn)*sqrt(2*E - vx0.^2); + + v = (vx0 + vy0); + KE = 0.5*(vx0.^2 + v.^2); + + x0 = 0.0; % 0.1 + y0 = -0.0; % -0.2 + + y0 = [x0 y0 vx0 vy0]; + + tspan = [1 3000]; % 500 + + options = odeset('RelTol',1e-9); + + [t,y] = ode45(@f1,tspan,y0,options); + + + siz = size(t); +% y1 = mod(y(1:siz,1) + pi,2*pi)-pi; % angle-1 +% y2 = mod(y(1:siz,2) + pi,2*pi)-pi; % angle-2 +% y3 = mod(y(1:siz,3) + pi,2*pi)-pi; % angular speed-1 +% y4 = mod(y(1:siz,4) + pi,2*pi)-pi; % antular speed-2 + y1 = y(1:siz,1); % angle-1 + y2 = y(1:siz,2); % angle-2 + y3 = y(1:siz,3); + y4 = y(1:siz,4); + + figure(1) + plot(t(1:siz),y1,t(1:siz),y2) + + figure(2) + plot(t(1:siz),y3,t(1:siz),y4) + xlabel('Time') + legend('speed') + + figure(3) + plot(y1,y2) + xlabel('x') + ylabel('y') + + KinEnergy = 0.5*y3.^2 + 0.5*(y3.^2 + y4.^2 + 2*y3.*y4.*cos(y2-y1)); + PotEnergy = - 2*cos(y1) - cos(y2); + Energy = KinEnergy + PotEnergy; + figure(4) + plot(t(1:siz),KinEnergy,t(1:siz),PotEnergy,t(1:siz),Energy) + title('Energy') + legend('Kin','Pot','E') + + % Power Spectrum + Pow = 0; + if Pow == 1 + yf = y(:,1); + F = fft(yf); + Pow = F.*conj(F); + + figure(4) + plot(Pow) + end + + trig = y2; % Trigger for crossing Poincare section: try y1 and switch xvar/yvar and vx/vy + + %First Return Map + Fst = 1; + if Fst == 1 + cnt = 0; + last = trig(1); + for loop = 2:siz + if (last < 0)&&(trig(loop) > 0) + cnt = cnt+1; + del1 = -trig(loop-1)./(trig(loop) - trig(loop-1)); + + vx(cnt) = y3(loop-1) + del1*(y3(loop)-y3(loop-1)); + vy(cnt) = y4(loop-1) + del1*(y4(loop)-y4(loop-1)); + xvar(cnt) = y1(loop-1) + del1*(y1(loop)-y1(loop-1)); + yvar(cnt) = y2(loop-1) + del1*(y2(loop)-y2(loop-1)); + + px(cnt) = 2*vx(cnt) + vy(cnt)*cos(yvar(cnt)-xvar(cnt)); + py(cnt) = vy(cnt) + vx(cnt)*cos(yvar(cnt)-xvar(cnt)); + + last = trig(loop); + else + last = trig(loop); + end + end + + figure(5) + plot(xvar,px,'o','MarkerSize',2,'Color',h(floor(0.9*mulnum*reploop)+1,:),'MarkerFaceColor',h(floor(0.9*mulnum*reploop)+1,:)) + xlabel('y') + ylabel('vy') + end + + pause(0.001) + +end % end reploop + +figure(5) +hold off +print -dtiff DoublePendPrint + + +% Model function + function dy = f1(t,y) + + A = y(4).^2.*sin(y(2)-y(1)); + B = -2*sin(y(1)); + C = y(3).^2.*sin(y(2)-y(1)).*cos(y(2)-y(1)); + D = sin(y(2)).*cos(y(2)-y(1)); + EE = 2 - (cos(y(2)-y(1))).^2; + + FF = y(4).^2.*sin(y(2)-y(1)).*cos(y(2)-y(1)); + G = -2*sin(y(1)).*cos(y(2)-y(1)); + H = 2*y(3).^2.*sin(y(2)-y(1)); + I = 2*sin(y(2)); + JJ = (cos(y(2)-y(1))).^2 - 2; + + dy = zeros(4,1); + dy(1) = y(3); + dy(2) = y(4); + dy(3) = (A + B + C + D)./EE; + dy(4) = (FF + G + H + I)./JJ; + + end % end f5 + + +end % end ltest + diff --git a/Heiles.m b/Heiles.m new file mode 100644 index 0000000..6343183 --- /dev/null +++ b/Heiles.m @@ -0,0 +1,161 @@ +function Heiles + +c(1) = 'r'; +c(2) = 'b'; +c(3) = 'g'; +c(4) = 'k'; +c(5) = 'c'; +c(6) = 'm'; +c(7) = 'y'; + +h = colormap(lines); + + +E = 0.1; % 1 0.1 + +epsE = 1; % 0.3425 1 + +prms = sqrt(E); +pmax = sqrt(2*E); + +% Potential Function +for xloop = 1:100 + x = -2 + 4*xloop/100; + for yloop = 1:100 + y = -2 + 4*yloop/100; + + V(yloop,xloop) = 0.5*x^2 + 0.5*y^2 + epsE*(x^2*y - 0.33333*y^3); + + end +end +figure(6) +colormap(gray) +imagesc(V) +hold on +contour(V,30,'LineColor','k') +hold off +caxis([-0.5 1]) +colorbar + +%keyboard + +h2 = figure(2); +h2.Position = [11 -31 560 420]; +dum = set(h2); + +h3 = figure(3); +h3.Position = [10 381 560 420]; +dum = set(h3); + + +h5 = figure(5); +h5.Position = [229 1 1008 804]; +dum = set(h5); +clf +hold on +repnum = 256; % 32 +mulnum = 64/repnum; +for reploop = 1:repnum + + clear yvar py + +% px = 2*(rand-0.499)*pmax; +% yp1 = 2*(rand-0.499)*sqrt(2*(E-px^2/2)); +% py = 2*(rand-0.499)*sqrt(2*(E-px^2/2-yp1^2/2)); +% xp1 = sign(rand-0.5)*real(sqrt(2*(E-px^2/2-yp1^2/2-py^2/2))); + + %Etest = px^2/2 + py^2/2 + xp1^2/2 + yp1^2/2 + + px = 2*(rand-0.499)*pmax; + py = sign(rand-0.5)*real(sqrt(2*(E-px^2/2))); + xp1 = 0; + yp1 = 0; + + y0 = [xp1 yp1 px py]; + + + + tspan = [1 2000]; % 500 + +% figure(1) + %options = odeset('OutputFcn',@odephas3); + options = odeset('RelTol',1e-9); + + [t,y] = ode45(@f1,tspan,y0,options); + + siz = size(t); + y1 = y(:,1); + y2 = y(:,2); + y3 = y(:,3); + y4 = y(:,4); + + figure(2) + plot(t(1:siz),y(1:siz,1)) + xlabel('Time') + legend('speed') + + figure(3) + plot(y(1:siz,1),y(1:siz,2)) + xlabel('x') + ylabel('y') + + + % Power Spectrum + Pow = 0; + if Pow == 1 + yf = y(:,1); + F = fft(yf); + Pow = F.*conj(F); + + figure(4) + plot(Pow) + end + + + %First Return Map + Fst = 1; + if Fst == 1 + cnt = 0; + last = y1(1); + for loop = 2:siz + if (last < 0)&(y1(loop) > 0) + cnt = cnt+1; + del1 = -y1(loop-1)./(y1(loop) - y1(loop-1)); + py(cnt) = y4(loop-1) + del1*(y4(loop)-y4(loop-1)); + yvar(cnt) = y2(loop-1) + del1*(y2(loop)-y2(loop-1)); + last = y1(loop); + else + last = y1(loop); + end + end + + figure(5) + plot(yvar,py,'o','MarkerSize',2,'Color',h(floor(0.9*mulnum*reploop)+1,:),'MarkerFaceColor',h(floor(0.9*mulnum*reploop)+1,:)) + xlabel('y') + ylabel('py') + end + + pause(0.001) + +end % end reploop + +figure(5) +hold off +print -dtiff HeilesPrint + + +% Model function + function dy = f1(t,y) + + B = 0.00; %0.002 + dy = zeros(4,1); + dy(1) = y(3); + dy(2) = y(4); + dy(3) = -y(1) - epsE*(2*y(1)*y(2) + B*y(3)); + dy(4) = -y(2) - epsE*(y(1)^2 - y(2)^2 + B*y(4)); + + end % end f5 + + +end % end ltest + diff --git a/NaK.m b/NaK.m new file mode 100644 index 0000000..cf5aaee --- /dev/null +++ b/NaK.m @@ -0,0 +1,89 @@ + +function NaK + +I = 50; % 10 +gL = 8; % 8 +EL = -78; % -80 -78 +gNa = 20; % 20 +ENa = 60; % 60 +gK = 10; % 10 +EK = -90; %-90 +C = 1; % 1 +Vn = -25; % -25 -45 +kn = 5; % 5 +Vm = -20; % -20 +km = 15; % 15 +tau = 1; % 1 0.152 + +%I = 12; EL = -78; Vn = -42.5; tau = 1; % pg. 117 supercritical Hopf bifurcation +%I = 4.53; EL = -80; Vn = -25; tau = 1; % pg. 113 Homoclinic +I = 5; EL = -80; Vn = -25; tau = 0.152; % pg. 110 Bistability Try I = 4 vs. I = 5 y0=[-20,0.4],[-70,0.1] +%I=45;gL=1;gNa=4;gK=4;Vm=-30;km=7;EL=-78;Vn=-45;tau=1; % pg. 176 subcritical Hopf + + +y0 = [-61,0.001]; +%y0 = [-20 0.4]; +%y0 = [10 0.81]; +%y0 = [-50 0.25]; + + + +tspan = [0 100]; + +figure(1) +options = odeset('OutputFcn',@odephas2); +[t,y] = ode45(@f5,tspan,y0,options); + +figure(2) +plot(t,y(:,1),t,y(:,2)) + + +figure(3) +plot(y(:,1),y(:,2)) +set(gcf, 'color', 'white') + +for loop = 1:100 + xp(loop) = -80 + 100*loop/100; + ninf = 1./(1 + exp((Vn - xp(loop))/kn)); + minf = 1./(1 + exp((Vm - xp(loop))/km)); + yp(loop) = ninf; + yyp(loop) = (I - gL*(xp(loop) - EL) - gNa*minf*(xp(loop) - ENa))/(gK*(xp(loop) - EK)); +end + +hold on +plot(xp,yp,'r',xp,yyp,'k') +hold off + +printfile = 0; + +if printfile == 1 + [sz,dum] = size(t); + ty1 = y(1:300,1)'; + ty2 = y(1:300,2)'; + tyV = y(1:sz,1); + Printfile3('Nullclines.txt',xp,yp,yyp); + Printfile2('Trajectory.txt',ty1,ty2); + Printfile2('TimeSeries.txt',t',tyV'); +end + + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + function yd = f5(t,y) + + ninf = 1./(1 + exp((Vn - y(1))/kn)); + minf = 1./(1 + exp((Vm - y(1))/km)); + + yp(1) = I/C - gL*(y(1) - EL)/C - gNa*minf*(y(1) - ENa)/C - gK*y(2).*(y(1) - EK)/C; + yp(2) = (ninf - y(2))./tau; + + + + yd = [yp(1);yp(2)]; + + end % end f5 + +end + diff --git a/PredPrey.m b/PredPrey.m new file mode 100644 index 0000000..c9a6e34 --- /dev/null +++ b/PredPrey.m @@ -0,0 +1,47 @@ +% PredPrey.m +% 9-3-11 +% + + +function PredPrey + +alpha = 3; +beta = 1; +gamma = 3; %-1 +delta = 1; +epsilon = 0; +zeta = 0; %+1 + +xinf = delta/epsilon; +yinf = alpha/beta; +eps = .1; + +options = odeset('RelTol',1e-6,'AbsTol',1e-7); + + +%y0 = [(xinf + eps) (yinf + 0.0)]; + +y0 = [2 2]; + +tspan = [0 200]; + +[t,y] = ode45(@f5,tspan,y0); + +figure(1) +plot(t,y(:,1),t,y(:,2)) + +figure(2) +plot(y(:,1),y(:,2)) +axis([0 10 0 10]) + + function yd = f5(t,y) + + yp(1) = y(1)*(alpha - beta*y(2) + epsilon*y(1)); + yp(2) = -y(2) *(gamma - delta*y(1) + zeta*y(2)); + + yd = [yp(1);yp(2)]; + + end % end f5 + +end + diff --git a/avgeigs.m b/avgeigs.m new file mode 100644 index 0000000..d6b2fc0 --- /dev/null +++ b/avgeigs.m @@ -0,0 +1,50 @@ + +N = 50; + +Nrep = 100; +eigvalER = zeros(1,N); +eigvalSF = zeros(1,N); +eigvalSW = zeros(1,N); +eigER = zeros(1,N); +eigSF = zeros(1,N); +eigSW = zeros(1,N); +for loop = 1:Nrep + + nodeER = makeER(N,0.08); + nodeSF = makeSF(N,2); + nodeSW = makeSW(N,2,0.1); + + [AER,degree,LapER] = adjacency(nodeER); + [VER,DER] = eig(LapER); + [vER,dER] = eig(AER); + [ASF,degree,LapSF] = adjacency(nodeSF); + [VSF,DSF] = eig(LapSF); + [sSF,dSF] = eig(ASF); + [ASW,degree,LapSW] = adjacency(nodeSW); + [VSW,DSW] = eig(LapSW); + [vSW,dSW] = eig(ASW); + for nloop = 1:N + eigvalER(nloop) = eigvalER(nloop) + DER(nloop,nloop)/Nrep; + eigvalSF(nloop) = eigvalSF(nloop) + DSF(nloop,nloop)/Nrep; + eigvalSW(nloop) = eigvalSW(nloop) + DSW(nloop,nloop)/Nrep; + eigER(nloop) = eigER(nloop) + dER(nloop,nloop)/Nrep; + eigSF(nloop) = eigSF(nloop) + dSF(nloop,nloop)/Nrep; + eigSW(nloop) = eigSW(nloop) + dSW(nloop,nloop)/Nrep; + end + +end + +x = 1:N; +figure(1) +plot(x,eigvalER,x,eigvalSF,x,eigvalSW) +legend('ER','SF','SW') +title('Graph Laplacian') +Printfile4('GLeig',x,eigvalER,eigvalSF,eigvalSW) + +figure(2) +plot(x,eigER,x,eigSF,x,eigSW) +legend('ER','SF','SW') +title('Adjacency') +Printfile4('Aeig',x,eigER,eigSF,eigSW) + + diff --git a/continued.m b/continued.m new file mode 100644 index 0000000..2636efa --- /dev/null +++ b/continued.m @@ -0,0 +1,20 @@ +% a = continued(r,N) +% continued fraction of rational number r with N terms up to N = 38 +% see also convergent.m + +function a = continued(r,N) + +rtmp = r; + +for loop = 1:N + + tmp = fix(rtmp); + if abs(rtmp - tmp) > 0 + rtmp = 1./(rtmp - tmp); + else + rtmp = 0; + end + + a(loop) = tmp; + +end \ No newline at end of file diff --git a/convergent.m b/convergent.m new file mode 100644 index 0000000..496edcd --- /dev/null +++ b/convergent.m @@ -0,0 +1,40 @@ +%function [p q] = convergent(r,N) +% Nth convergent of real number r (up to N = 38) +% p and q are relatively prime +% +% see also continued.m + + +function [p q] = convergent(r,Ntmp) + +N = Ntmp + 1; +a = continued(r,N); + +num = 1; +denom = a(N); + +flag = 1; +while flag == 1 + if denom == 0 + N = N-1; + denom = a(N); + else + flag = 0; + end +end + +for loop = 1:N-1 + + if(loop<(N-1)) + num = a(N-loop)*denom + num; + tmp = denom; + denom = num; + num = tmp; + else + num = a(N-loop)*denom + num; + end +end + +p = num; +q = denom; + diff --git a/coupledrossler.m b/coupledrossler.m new file mode 100644 index 0000000..27a2a4a --- /dev/null +++ b/coupledrossler.m @@ -0,0 +1,82 @@ +% function rossler + +function y = coupledrossler + +a = 0.25; % 0.15 (defined phase) 0.3 (chaotic phase) +b = 0.4; % 0.4 +c = 8.5; % 6 8.0 +d = 0.97; +g = 0.1; + +y0 = [1 1 0 0 1 0]; +tspan = [1 30]; +[t,y] = ode45(@f5,tspan,y0); +[sz dum] = size(y); + +y0 =[y(sz,1) y(sz,2) y(sz,3) y(sz,4) y(sz,5) y(sz,6)]; +tspan = [1 200]; +options = odeset('OutputFcn',@odephas3); +[t,y] = ode45(@f5,tspan,y0,options); + +figure(2) +plot(t,y(:,1),t,y(:,2),t,y(:,3)) + +figure(3) +plot(y(:,1),y(:,6)) +title('x1-z2') + +figure(4) +plot(y(:,2),y(:,6)) +title('y1-z2') + +figure(5) +plot(y(:,1),y(:,5)) +title('x1-y2') + +theta1 = atan2(y(:,2),y(:,1)); +theta2 = atan2(y(:,5),y(:,4)); + +theta12 = atan2(y(:,2),y(:,4)); +theta21 = atan2(y(:,5),y(:,1)); + +figure(6) +plot(theta1,theta2) +title('theta1 vs. theta2') + +figure(7) +plot(theta21-theta12) +title('theta21-theta12') + +phse = theta1 - theta2; +figure(8) +[temp bin] = maskbilevel(phse,-pi,pi,0,0); +plot(temp) +title('theta1-theta2') + +Metric1 = sqrt(sum((temp.*bin).^2))/sqrt(sum(bin.^2)) + + + +d = y(:,1); +e = y(:,2); +f = y(:,3); +Printfile4('lout',t',d',e',f'); + + + function yd = f5(t,y) + + yp(1) = -d*y(2) - y(3) + g*(y(4) - y(1)); + yp(2) = d*y(1) + a* y(2) + g*(y(5) - y(2)); + yp(3) = b + y(3)*(y(1) - c) + g*(y(6) - y(3)); + + yp(4) = -1*d*y(5) - y(6) + g*(y(1) - y(4)); + yp(5) = 1*d*y(4) + a* y(5) + g*(y(2) - y(5)); + yp(6) = b + y(6)*(y(4) - c) + g*(y(3) - y(6)); + + yd = [yp(1);yp(2);yp(3);yp(4);yp(5);yp(6)]; + + end % end f5 + + +end % end ltest + diff --git a/diffusiondriver.m b/diffusiondriver.m new file mode 100644 index 0000000..fab31c1 --- /dev/null +++ b/diffusiondriver.m @@ -0,0 +1,158 @@ +% diffusiondriver.m +% 3-24-12 + +clear + +close all + +N = 128; +p = 0.05; +m = 4; +beta = 0.01; + + +node = makeER(N,0.06); +%node = makeSF(N,m); +%node = makeSW(N,m,0.1); + + +% A = adjacency(node); +% A = ham2adj(N); +% node = adj2node(A); + +[N,e,avgdegree,maxdegree,mindegree,numclus,meanclus,Lmax,L2,LmaxL2] = clusterstats(node); + +disp(' ') +displine('Number of nodes = ',N) +disp(strcat('Number of edges = ',num2str(e))) +disp(strcat('Mean degree = ',num2str(avgdegree))) +displine('Maximum degree = ',maxdegree) +disp(strcat('Number of clusters = ',num2str(numclus))) +disp(strcat('mean cluster coefficient = ',num2str(meanclus))) +disp(' ') +disp(strcat('Lmax = ',num2str(Lmax))) +disp(strcat('L2 = ',num2str(L2))) +disp(strcat('Lmax/L2 = ',num2str(LmaxL2))) +disp(' ') + + [A,degree,Lap] = adjacency(node); + + [V,D] = eig(Lap); + + for loop = 1:N + eigval(loop) = D(loop,loop); + end + + figure(1) + plot(eigval) + title('Eigenvalues') + + % initial values + c = zeros(N,1); + c(1) = 1; + + % eigvec decomposition + + for eigloop = 1:N + Vtemp = V(:,eigloop); + v(eigloop) = sum(c.*Vtemp); + end + + % time loop + Ntime = 100; + for timeloop = 1:Ntime % 200 + + for nodeloop = 1:N + + temp = 0; + for eigloop = 1:N + + temp = temp + V(nodeloop,eigloop)*v(eigloop)*exp(-eigval(eigloop)*beta*(timeloop-1)); + + end % end eigloop + + concentration(timeloop,nodeloop) = temp; + + end % endnodeloop + + + end % end timeloop + + figure(2) + imagesc(real(log(concentration))) + colormap(jet) + colorbar + caxis([-10 0]) + title('Log Concentrations vs. time') + xlabel('Node Number') + + figure(3) + plot(concentration(100,:)) + title('Ending Concentrations') + xlabel('Node Number') + + + x = 0:Ntime-1; + h = colormap(jet); + figure(4) + for nodeloop = 1:N + rn = round(rand*63 + 1); + y = concentration(:,nodeloop)+0.001; + semilogy(x,y,'Color',h(rn,:)) + hold on + end + hold off +title('Concentrations vs. time') + + + x = 0:Ntime-1; + h = colormap(jet); + figure(5) + for nodeloop = 1:10 + rn = round(rand*63 + 1); + y = concentration(:,nodeloop*10)+0.001; + %semilogy(x,y,'Color',h(rn,:),'LineWidth',1.1) + semilogy(x,y,'k','LineWidth',1.2) + hold on + end + hold off + set(gcf,'Color','White') + title('Selected Nodes: Continuous time') + + + % Now try the discrete-time-map approach + + c0 = c; + dt = 1; % 5 + M = eye(N,N) - beta*Lap*dt; + + for timeloop = 1:200 %40 + + c = (M^timeloop)*c0; + + Con(timeloop,:) = c'; + end + + x = 0:1:199; % 0:5:199 + h = colormap(jet); + figure(6) + for nodeloop = 1:10 + rn = round(rand*63 + 1); + y = Con(:,nodeloop*10)+0.001; + %semilogy(x,y,'Color',h(rn,:),'LineWidth',1.1) + semilogy(x,y,'k','LineWidth',1.2) + hold on + end + hold off + set(gcf,'Color','White') + title('Selected Nodes: Discrete time') + + + + + + + + + + \ No newline at end of file diff --git a/flowsimple.m b/flowsimple.m new file mode 100644 index 0000000..5e9a8fa --- /dev/null +++ b/flowsimple.m @@ -0,0 +1,108 @@ +function flowsimple + +xrange = [-5 5]; +yrange = [-5 5]; +% xrange = [0 4]; +% yrange = [0 4]; +rngx = xrange(2) - xrange(1); +rngy = yrange(2) - yrange(1); + +[X,Y] = meshgrid(xrange(1):0.1:xrange(2), yrange(1):0.1:yrange(2)); + +smallarrow = 0; +[x,y] = f5(X,Y); + +clf +figure(1) +for xloop = 1:10 + xs = xrange(1) +xloop*rngx/10; + for yloop = 1:10 + ys = yrange(1) +yloop*rngy/10; + + streamline(X,Y,x,y,xs,ys) + + end +end +hold on +[XQ,YQ] = meshgrid(xrange(1):1:xrange(2),yrange(1):1:yrange(2)); +smallarrow = 0; +[xq,yq] = f5(XQ,YQ); +quiver(XQ,YQ,xq,yq,.2,'r','filled') +hold off + +axis([xrange(1) xrange(2) yrange(1) yrange(2)]) +set(gcf,'Color','White') + + + function [x,y] = f5(X,Y) + +% xtmp = X*(1-X^2) - Y; +% ytmp = X; + +% xtmp = Y; % Double well +% ytmp = -0.5*Y + X - X.^3; + +% xtmp = X*(1-2*X+Y); +% ytmp = Y*(1+X-Y); + + +% xtmp = X.*(1-Y.^3); % testing dsge +% ytmp = Y.*(2.5-X); + +% xtmp = X.*Y + 1 - 0.5*X; +% ytmp = Y.*(2 - X); + + xtmp = X.*(3-Y); % Rabbit sheep + ytmp = Y.*(2-X); + +% xtmp = Y; +% ytmp = -X - 0.1*exp(-abs(X)/1)*X.^2; + + +% c = 1; % Medio +% xtmp = -Y + X.*(c - Y.^2); +% ytmp = X + 1; + +% xtmp = X.*(1-2*X+Y); +% ytmp = Y.*(0.5 + 0.5*X - Y); + +% xtmp = Y+0.75*X; +% ytmp = X.^3-Y-X; + + +% xtmp =X.*(2-X-Y); +% ytmp = X-Y; + + +% xtmp = X.*(X-Y); +% ytmp = Y.*(2*X-Y); + + +% a = 2; b = 2; % Inward spiral +% xtmp = Y - a*X; +% ytmp = -X - b*Y; + +% xtmp = -X.^2-Y; +% ytmp = X; + +% xtmp = X.^2 - Y; +% ytmp = X - Y; + + %x = X.*(3-X-2*Y); + %y = Y.*(2-X-Y); + + nrm = sqrt(xtmp.^2 + ytmp.^2); + + if smallarrow == 1 + x = xtmp./nrm; + y = ytmp./nrm; + else + x = xtmp; + y = ytmp; + end + + end % end f5 + + +end % end flowsimple + diff --git a/logistic.m b/logistic.m new file mode 100644 index 0000000..2360df4 --- /dev/null +++ b/logistic.m @@ -0,0 +1,40 @@ +% logistic.m +% Logistic map + +clear + +xold = 0.231; +maxg = 900; % 4000 for good coverage (pause off) +ming = 2.9; +figure(1); +clf +axis([2.7 4 0 1]); +hold on +for gainloop = 1:maxg + + g=(4.00-ming)*gainloop/maxg + ming; + + for sloop = 1:100 %settle-down loop + xnew = g*xold*(1-xold); + xold = xnew; + end % end settle-down loop + + rep = 64; + ind = 0; + for rloop = 1:rep + xnew = g*xold*(1-xold); + xold = xnew; + + ind = ind + 1; + x(ind) = g + (4.00-ming)*(rloop-1)/(maxg*rep); + y(ind) = xnew; + + end % end rloop + + plot(x,y,'o','MarkerSize',2); + pause(0.01) % Take pause off to watch real-time, but put maxg to 300 + + +end % end gainloop +hold off + diff --git a/lorenz.m b/lorenz.m new file mode 100644 index 0000000..ff35dd0 --- /dev/null +++ b/lorenz.m @@ -0,0 +1,110 @@ +function lorenz + +a = 10; % 10 +b = 50; % 50 +c = 8/3; % 8/3 + +y0 = [20 1 1]; +%y0 = [12.8 10.955 53.455]; +%y0 = [-4.3426 -9.3985 21.389]; + +tspan = [0 20]; +[t,yy] = ode45(@f5,tspan,y0); +[sz,~] = size(yy); + +y0 = yy(sz,:); +tspan = [0 100]; + +options = odeset('OutputFcn',@odephas3); +% options = odeset('OutputFcn',@odeplot); +options = odeset('RelTol',1e-5); + +[t,y] = ode45(@f5,tspan,y0,options); + +figure(2) +colormap(jet) +plot(t,y(:,1),'k',t,y(:,2),'r',t,y(:,3),'b') +legend('x','y','z') + +figure(20) +colormap(jet) +plot3(y(:,1),y(:,2),y(:,3),'k') +set(gcf,'Color','white') +grid on + + +figure(3) +plot(y(:,1),y(:,2)) +xlabel('x-value') +ylabel('y-value') + +u = sqrt(y(:,1).^2 + y(:,2).^2); +z = y(:,3); + +figure(4) +plot(u,z) +xlabel('u-value') +ylabel('z-value') + +xp = u - 18; +yp = z - 52; + +theta = atan2(yp,xp); +signal = sin(theta); + +figure(5) +plot(signal,'LineWidth',2) + + +PowS = 1; +% Power Spectrum +if PowS == 1 % Turn this off by setting PowS = 0 + [st dum] = size(t); + tshort = t(400:st,1); + %yshort = y(400:st,1); + yshort = u(400:st,1); + [st dum] = size(tshort); + dt = (max(tshort)-min(tshort))/st; + df = 1/st/dt; + fmax = 1/dt; + freq = df:df:fmax; + + yf = yshort(:,1); + F = fft(yf); + Pow = F.*conj(F); + + mid = round(st/4); + figure(6) + semilogx(freq(2:mid),Pow(2:mid)) + title('Power Spectrum') + xlabel('Frequency (Hz)') + + + [Pmax I] = max(Pow); + f0 = freq(I) + T0 = 1/f0 +% T0nat = 2*pi/sqrt(w02) +end + + + + +d = y(:,1); +e = y(:,2); +f = y(:,3); +Printfile4('lout',t',d',e',f') + + + function yd = f5(t,y) + + yp(1) = a*(y(2) - y(1)); + yp(2) = y(1)*(b - y(3)) - y(2); + yp(3) = y(1)*y(2) - c*y(3); + + yd = [yp(1);yp(2);yp(3)]; + + end % end f5 + + +end % end ltest + diff --git a/lotkavoltstream.m b/lotkavoltstream.m new file mode 100644 index 0000000..f90e66a --- /dev/null +++ b/lotkavoltstream.m @@ -0,0 +1,97 @@ +% function lotkavoltstream +% + + +function lotkavoltstream + +alpha = 1; +beta = 1.5; +gamma = 1.; +delta = 1.25; + +figure(1) +N = 20; +for loop = 2:N + + xstar = gamma/delta; + ystar = alpha/beta; + + valx = xstar*loop/N; + valy = ystar*loop/N; + y0 = [valx valy]; + tspan = [0 10]; + [t,y] = ode45(@F5,tspan,y0); + [sz dum] = size(t); + y0 = [y(sz,1) y(sz,2)]; + clear t y + tspan = [10 20]; + options = odeset('RelTol',1e-7); + [t,y] = ode45(@F5,tspan,y0,options); + + plot(y(:,1),y(:,2),'LineWidth',1.5) + mxlim = 3*max(xstar,ystar); + axis([0 mxlim 0 mxlim]) + + hold on + +end +hold off +set(gcf,'Color','white') + + +figure(2) +N = 3; +for loop = 2:N + + xstar = gamma/delta; + ystar = alpha/beta; + + valx = xstar*(loop-1)/N; + valy = ystar*(loop-1)/N; + y0 = [valx valy]; + tspan = [0 10]; + [t,y] = ode45(@F5,tspan,y0); + [sz dum] = size(t); + y0 = [y(sz,1) y(sz,2)]; + clear t y + tspan = [10 40]; + options = odeset('RelTol',1e-7); + [t,y] = ode45(@F5,tspan,y0,options); + + plot(t,y(:,1),t,y(:,2),'--','LineWidth',1.5) + mxlim = 3*max(xstar,ystar); + + hold on + +end +hold off +set(gcf,'Color','white') + + + + + +%keyboard + + function [x,y] = f5(X,Y) + + + x = X.*(alpha - beta*Y); + y = -Y*(gamma - delta*X); + + + end % end f5 + + function yd = F5(t,y) + + yp(1) = y(1).*(alpha - beta*y(2)); + yp(2) = -y(2).*(gamma - delta*y(1)); + + + yd = [yp(1);yp(2)]; + + end % end f5 + + +end % end flowsimple + diff --git a/quasiHam.m b/quasiHam.m new file mode 100644 index 0000000..cc5e622 --- /dev/null +++ b/quasiHam.m @@ -0,0 +1,145 @@ +% Quasispecies simulation. Includes mutation with Hamming distance. + +function quasiHam + +%close all +h = colormap(lines); + +randpop = 0; + +B = 8; +N = 2^B; % size of mutation space (32) +ep = 0.1; % average mutation rate: 0.1 to 0.01 typical +lam = 1; % gain of Hamming distance + +if randpop == 1 + x0temp = rand(1,N); % Initial population + sx = sum(x0temp); + x0 = x0temp/sx; +else + x0 = zeros(1,N); + x0(1) = 0.667; x0(2) = 0.333; +end + +Pop0 = sum(x0); + +for yloop = 1:N + for xloop = 1:N + H(yloop,xloop) = hamming(yloop,xloop); + end +end + +%keyboard + + +Qtemp = 1./(1+ep*H); %Mutation matrix +Qsum = sum(Qtemp,2); + +% Normalize along species +for yloop = 1:N + for xloop = 1:N + Q(yloop,xloop) = Qtemp(yloop,xloop)/Qsum(yloop); + end +end + +xp = 1:N; +alpha = 64; +ftemp = exp(-lam*H(alpha,:)); % Fitness landscape +sf = sum(ftemp); +f = ftemp/sf; + +% Transition matrix +for yloop = 1:N + for xloop = 1:N + W(yloop,xloop) = f(yloop)*Q(yloop,xloop); + end +end + +tspan = [0 500]; + + +[t,x] = ode45(@quasispec,tspan,x0); + +Pop0 +[sz,dum] = size(t); +Popend = sum(x(sz,:)) + +phistar = sum(f.*x(sz,:)) + + +figure(1) +plot(f,'o-') +hold on +figure(1) +plot(x(sz,:),'o-r') +hold off +legend('fitness','population') + +figure(2) +for loop = 1:N + semilogx(t,x(:,loop),'Color',h(1+round(loop*63/N),:)) + hold on +end +hold off + +figure(3) +for loop = 1:N + plot(t,x(:,loop),'Color',h(1 + round(loop*63/N),:),'LineWidth',1.25) + hold on +end +hold off +set(gcf,'Color','white') +xlabel('Time','FontSize',14) +ylabel('Population','FontSize',14) +hh = gca; +set(hh,'FontSize',14) + +figure(4) +for loop = 1:N + loglog(t,x(:,loop),'Color',h(1+round(loop*63/N),:)) + hold on +end +hold off + +figure(5) +for loop = 1:N + semilogy(t,x(:,loop),'Color',h(1+round(loop*63/N),:),'LineWidth',1.25) + hold on +end +hold off +set(gcf,'Color','white') +xlabel('Time','FontSize',14) +ylabel('Population','FontSize',14) +hh = gca; +set(hh,'FontSize',14) + +%keyboard + + + +% Eigenvalues + +[V,D] = eig(W); + +max(D(:,1)) + +figure(6) +plot(V(:,2)) +title('eigenvalues') + + +disp(' ') + + + function yd = quasispec(t,y) + + phi = sum(f'.*y); % Average fitness of population + + yd = W*y - phi*y; + + + end % end quasispec + +end % end quasi + + diff --git a/quasiSpec.m b/quasiSpec.m new file mode 100644 index 0000000..f1a511e --- /dev/null +++ b/quasiSpec.m @@ -0,0 +1,223 @@ +% Quasispecies simulation. Includes mutation with Hamming distance. + +function quasiSpec + +%close all +global W +h = colormap(lines); + +randpop = 0; % 0) = spike population; 1) = random population +mutype = 0; % 0) = Hamming; 1) = rand +fitype = 0; % 0) = Hamming; 1) = 2-peak; 2) = rand+gauss; 3) freq-dep + +B = 7; +N = 2^B; % size of mutation space (64) +lam = 1; % Hamming fitness only +gamma = 1; % freq-dep fitness (payoff matrix only) +relran = 0.025; % relative random contrib to fitness +time_expand = 50; + +ep = 0.08; % average mutation rate: 0.1 to 0.01 typical (0.4835) (0.0290) + + +%%%%% Set original population +if randpop == 1 + rng(0); + x0temp = rand(1,N); % Initial population + sx = sum(x0temp); + x0 = x0temp/sx; +else + x0 = zeros(1,N); + x0(1) = 0.667; x0(2) = 0.333; +end + +Pop0 = sum(x0); + + +%%%%%% Set Hamming distance +for yloop = 1:N + for xloop = 1:N + H(yloop,xloop) = hamming(yloop-1,xloop-1); + end +end + + + +%%%%%%% Set Mutation matrix +if mutype == 0 + Qtemp = 1./(1+H/ep); %Mutation matrix on Hamming + %Qtemp = exp(-H/(ep*50)); + Qsum = sum(Qtemp,2); + + % Normalize mutation among species + for yloop = 1:N + for xloop = 1:N + Q(yloop,xloop) = Qtemp(yloop,xloop)/Qsum(xloop); + end + end + +end +if mutype == 1 + rng(0); + S = stochasticmatrix(N); + Stemp = S - diag(diag(S)); + Qtemp = ep*Stemp; + sm = sum(Qtemp,2)'; + Q = Qtemp + diag(ones(1,N) - sm); +end + + +%keyboard + + +%%%%%%% Set fitness landscape + +if fitype == 0 % Hamming + x = 1:N; + alpha = 84; + ftemp = exp(-lam*H(alpha,:)); % Fitness landscape + sf = sum(ftemp); + f = ftemp/sf; +end +if fitype == 1 % double peak and rand + rng(1); + f = rand(1,N); + x = 1:N; + delg = 20; + sig1 = 1; + sig2 = 4; + g1 = gaussprob(x,(N/2 - delg),sig1); + g2 = 3*gaussprob(x,(N/2 + delg),sig2); + ftemp = relran*f + g1 + g2; + f = ftemp/sum(ftemp); +end +if fitype == 2 % rand + Gauss + rng(0); + f = rand(1,N); + x = 1:N; + ftemp = relran*f + gauss((x-N/2)/2); % Fitness landscape + f = ftemp/sum(ftemp); +end +if fitype == 3 % frequency-dependent Hamming + avgdis = mean(mean(H)); + %payoff = exp(-gamma*(H - avgdis)); % payoff matrix + %payoff = H.^2; + %payoff = ones(size(H)); + payoff = exp(-gamma*H); +end + +%keyboard + +% Run time evolution +tspan = [0 1000]; +[t,x] = ode45(@quasispec,tspan,x0); + +Pop0 +[sz,dum] = size(t); +Popend = sum(x(sz,:)) + +phistar = sum(f.*x(sz,:)) % final average fitness + + +figure(1) +plot(f,'-') +hold on +figure(1) +plot(x(sz,:),'r') +hold off +legend('fitness','population') + + +figure(2) +for loop = 1:N + semilogx(t,x(:,loop),'Color',h(round(loop*64/N),:),'LineWidth',1.25) + hold on +end +hold off +set(gcf,'Color','white') +xlabel('Time','FontSize',14) +ylabel('Population','FontSize',14) +hh = gca; +set(hh,'FontSize',14) + + +figure(3) +for loop = 1:N + plot(t,x(:,loop),'Color',h(round(loop*64/N),:)) + hold on +end +hold off + +figure(4) +for loop = 1:N + loglog(t,x(:,loop),'Color',h(round(loop*64/N),:)) + hold on +end +hold off + +figure(5) +for loop = 1:N + semilogy(t,x(:,loop),'Color',h(round(loop*64/N),:)) + hold on +end +hold off + + + +% Eigenvalues + +[V,D] = eig(W); + +max(D(:,1)) + +figure(6) +%semilogy(abs(V(:,1))) +plot((V(:,1))) +title('Eigenvector') + +disp(' ') + + +if fitype == 1 + xlo = N/2 - delg - 2*sig1; + xhi = N/2 - delg + 2*sig1; + fit44 = sum(f(xlo:xhi)) + pop44 = sum(x(sz,xlo:xhi))/Popend + xlo = N/2 + delg - 2*sig2; + xhi = N/2 + delg + 2*sig2; + fit84 = sum(f(xlo:xhi)) + pop84 = sum(x(sz,xlo:xhi))/Popend + +end + + + + function yd = quasispec(~,y) + + if fitype == 3 % frequency-dependent Hamming + for loop = 1:N + ftemp(loop) = sum(payoff(:,loop).*y); + end + f = time_expand*ftemp/sum(ftemp); + end + + + % Transition matrix + for yloop = 1:N + for xloop = 1:N + W(yloop,xloop) = f(yloop)*Q(yloop,xloop); + end + end + + phi = sum(f'.*y); % Average fitness of population + + yd = W*y - phi*y; + + + end % end quasispec + + + + +end % end quasiSpec + diff --git a/randsgn.m b/randsgn.m new file mode 100644 index 0000000..a3155d2 --- /dev/null +++ b/randsgn.m @@ -0,0 +1,10 @@ +% function y = randsgn(N) + + +function y = randsgn(N) + +%x = randbin(N,0.5); + +ytmp = 2*(0.5-rand(1,N)); +y = sign(ytmp); + diff --git a/repeq.m b/repeq.m new file mode 100644 index 0000000..f4f30aa --- /dev/null +++ b/repeq.m @@ -0,0 +1,194 @@ +% replicator equation +% 3-10-12 + +function repeq + +N = 7; +asymm = 4; % 1 = zero diag (replicator eqn) 2 = zero trace (autocatylitic model) 3 = winner take all 4 = hypercycle +phi0 = 0.01; % average fitness (positive number) damps oscillations +fac = 1; + +displine('N = ',N) +displine('asymm = ',asymm) +displine('phi0 = ',phi0) + +%h = newcolormap('fluorodark'); +h = colormap(jet); +hg = [1 1 1] - colormap(gray); + +tempx = rand(1,N); +x0 = tempx/sum(tempx);% initial populations + +node = makeER(N,1); +%node = makeSW(N,N/32,0.25); +%node = makeSF(N,N/32); +Adj = adjacency(node); +numclus = clusternum(node); +[N,e,n] = clusterstats(node); +displine('number of clusters = ',numclus) +displine('average degree = ',n) + +if asymm == 1 + % Asymmetric payoff matrix with zero diagonals (replicator model) + A = zeros(N,N); + for yloop = 1:N + for xloop = yloop + 1:N + A(yloop,xloop) = 2*(0.5 - rand); + A(xloop,yloop) = -A(yloop,xloop); + end + end + +elseif asymm == 2 + % Asymmetric payoff matrix with zero trace (autocatylitic model) + Atemp = zeros(N,N); + for yloop = 1:N + for xloop = yloop+1:N + Atemp(yloop,xloop) = 2*(0.5 - rand); + Atemp(xloop,yloop) = -Atemp(yloop,xloop); + end + Atemp(yloop,yloop) = 2*(0.5 - rand); + end + tr = trace(Atemp); + A = Atemp; + + for yloop = 1:N + A(yloop,yloop) = Atemp(yloop,yloop) - tr/N; + end + +elseif asymm == 3 + % Alternating positive/negative payoff matrix with zero diagonals (winner take all model) + A = zeros(N,N); + for yloop = 1:N-1 + xloop = yloop + 1; + A(yloop,xloop) = odd(xloop)*abs(randn)-even(xloop)*abs(randn); + A(xloop,yloop) = -A(yloop,xloop); + end + +elseif asymm == 4 + % hypercycle + A = zeros(N,N); + for yloop = 2:N + A(yloop,yloop-1) = 0.1*N + rand; + end + A(1,N) = 0.1*N + rand; + + +end % end if asymm + +%keyboard + +E = eye(N,N); +A = fac*(A.*Adj + E.*A); + +%keyboard + +figure(1) +imagesc(A); +colormap(h) +colorbar +title('Payoff Matrix') + +tspan = [1 1600]; +[t,x] = ode45(@f5,tspan,x0); + +[sy,sx] = size(x); + +% fitness +xend = x(sy,:)'; +fit = A*xend; + +avgfit = fit'*xend; + +pop = sum(x,2); +av = mean(mean(x)); +displine('population = ',mean(pop)) +displine('av pop =', mean(mean(x))); + +% How many non-zero? +cnt = 0; +for loop = 1:N + if mean(x(sy-50:sy,loop),1) > av/20 + cnt = cnt + 1; + end +end +displine('cnt = ',cnt) +displine('avg fitness = ',avgfit) +xend +fit +A + +% Average max slope? +% for loop = 1:N +% mx = max(x(200:sy,loop)); +% deriv(:,loop) = diff(x(200:sy,loop)); +% freq(loop) = mean(abs(deriv(:,loop)))/mx; +% end +% mnfreq = mean(freq); +% displine('mean frequency = ',mnfreq) + +figure(2) +for loop = 1:N + %plot(t,real(log(x(:,loop))),'Color',h(round(loop*64/N),:),'LineWidth',1.25) + %axis([0 1600 -10 0]) + plot(t,x(:,loop),'Color',h(round(loop*64/N),:),'LineWidth',1.25) + hold on +end +hold off + +figure(3) +for loop = 1:N + plot(t,real(log(x(:,loop))),'Color',h(round(loop*64/N),:),'LineWidth',1.25) + axis([0 1600 -10 0]) + %plot(t,x(:,loop),'Color',h(round(loop*64/N),:),'LineWidth',1.25) + hold on +end +hold off + +% figure(3) +% plot(x(:,1),x(:,N)) +% +% figure(4) +% plot(t,x(:,1),t,(x(:,round(N/2))),t,x(:,N)) + +%keyboard + + + function yd = f5(t,y) + + for iloop = 1:N + ftemp = 0; + for jloop = 1:N + ftemp = ftemp + A(iloop,jloop)*y(jloop); + end % end jloop + f(iloop) = ftemp; + end % end iloop + + phitemp = phi0; % Can adjust this from 0 to 1 to stabilize (but Nth population is no longer independent) + for loop = 1:N + phitemp = phitemp + f(loop)*y(loop); + end + phi = phitemp; + + for loop = 1:N-1 + yd(loop) = y(loop)*(f(loop) - phi); + end + + if abs(phi0) < 0.01 % average fitness maintained at zero + yd(N) = y(N)*(f(N)-phi); + else % non-zero average fitness + ydtemp = 0; + for loop = 1:N-1 + ydtemp = ydtemp - yd(loop); + end + yd(N) = ydtemp; + end + + yd = yd'; + + + end % end f5 + + +end % end repeq + + diff --git a/repmut.m b/repmut.m new file mode 100644 index 0000000..f8d57fd --- /dev/null +++ b/repmut.m @@ -0,0 +1,164 @@ +%function repmut +% 12/08/14 + +function repmut + +clear +format compact + + +f = 1; % this is a dummy init to make f global +N = 64; % 64 +p = 1/sqrt(N); +time_expand = N; + +mutype = 0; % 0 = Hamming 1 = rand +pay = 0; % 0 = Hamming 1 = 1/sqrt(N) +ep = .1; % average mutation rate: 0.1 to 0.01 typical (0.4835) + +%%%%% Set original population +x0temp = rand(1,N); % Initial population +sx = sum(x0temp); +y0 = x0temp/sx; +Pop0 = sum(y0); + + +%%%%% Set Adjacency + +node = makeglobal(N); +%node = makeER(N,0.1); +%node = makeSF(N,4); +%node = makeSW(N,4,0.5); +[Adj,degree,Lap] = adjacency(node); + + + +%%%%%% Set Hamming distance +for yloop = 1:N + for xloop = 1:N + H(yloop,xloop) = hamming(yloop-1,xloop-1); + end +end + + + +%%%%%%% Set Mutation matrix +if mutype == 0 + Qtemp = 1./(1+H/ep); %Mutation matrix on Hamming + %Qtemp = exp(-H/(ep*50)); + Qsum = sum(Qtemp,2); + + % Normalize mutation among species + for yloop = 1:N + for xloop = 1:N + Q(yloop,xloop) = Qtemp(yloop,xloop)/Qsum(xloop); + end + end + +elseif mutype == 1 + S = stochasticmatrix(N); + Stemp = S - diag(diag(S)); + Qtemp = ep*Stemp; + sm = sum(Qtemp,2)'; + Q = Qtemp + diag(ones(1,N) - sm); +end + + +figure(1) +imagesc(Q) +title('Mutation Matrix') +colormap(jet) + +%%%%%%% Set payoff matrix +if pay == 1 + payoff = zeros(N,N); + for yloop = 1:N + payoff(yloop,yloop) = 1; + for xloop = yloop + 1:N + payoff(yloop,xloop) = p; + payoff(xloop,yloop) = p; + %payoff(yloop,xloop) = p*2*(0.5 - randbin(1,0.5)); + %payoff(xloop,yloop) = payoff(yloop,xloop); + %payoff(xloop,yloop) = -payoff(yloop,xloop); + end + end +elseif pay == 0 + payoff = exp(-1*H); +end + + +figure(2) +imagesc(payoff) +title('Payoff Matrix') +colormap(jet) + + +% Run time evolution +tspan = [0 1000]; +[t,x] = ode45(@quasispec,tspan,y0); + +Pop0 +[sz,dum] = size(t); +Popend = sum(x(sz,:)) + +phistar = sum(f.*x(sz,:)) % final average fitness + + + +figure(3) +clf +h = colormap(lines); +for loop = 1:N + plot(t,x(:,loop),'Color',h(round(loop*64/N),:)) + hold on +end +hold off + +figure(4) +clf +for loop = 1:N + semilogx(t,x(:,loop),'Color',h(round(loop*64/N),:)) + hold on +end +hold off + +figure(5) +clf +for loop = 1:N + semilogy(t,x(:,loop),'Color',h(round(loop*64/N),:)) + hold on +end +hold off + +%keyboard + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + function yd = quasispec(~,y) + + for loop = 1:N + ftemp(loop) = sum(payoff(:,loop).*y); + end + f = time_expand*ftemp/sum(ftemp); + + + % Transition matrix + for yloop = 1:N + for xloop = 1:N + W(yloop,xloop) = f(yloop)*(Adj(yloop,xloop)*Q(yloop,xloop)); + end + end + + phi = sum(f'.*y); % Average fitness of population + + yd = W*y - phi*y; + + + end % end quasispec + +end + + diff --git a/rksimple.m b/rksimple.m new file mode 100644 index 0000000..e3b9115 --- /dev/null +++ b/rksimple.m @@ -0,0 +1,53 @@ + +function rksimple + +clear + + +y0 = [1 0]; +%y0 = [4,4]; +tspan = [0 400]; + +[t,y] = ode45(@f5,tspan,y0); + +figure(1) +%plot(t,y(:,1),t,y(:,2),t,y(:,3)) +plot(t,y(:,1),t,y(:,2)); +legend('1','2','3') + +figure(2) +plot(y(:,1),y(:,2)) + + + function yd = f5(t,y) + + w = 1; + + yp(1) = y(2); + yp(2) = -y(1) + 0.5*y(2)*(1-y(1).^2); + + %mu = 1; beta = 1; w02 = 1; + %yp(1) = y(2); + %yp(2) = mu*y(2)*(1-beta*y(1)^2)-w02*y(1); + +% yp(1) = -y(2) - y(1).^3 + 5*y(1); +% yp(2) = -y(2) + 4*(y(1)-0.0); + +% mu = 1.5; beta = 1; +% yp(2) = -mu*y(2)+beta*y(1)*y(2); %SIR +% yp(1) = -beta*y(1)*y(2); + +% yp(1) = w*(y(2)-y(3)); +% yp(2) = w*(y(3)-y(1)); +% yp(3) = w*(y(1)-y(2)); + +% yp(1) = y(1)*(1 - 2*y(1) + y(2)); +% yp(2) = y(2)*(1 + y(1) - y(2)); + + %yd = [yp(1);yp(2);yp(3)]; + yd = [yp(1);yp(2)]; + + end % end f5 + +end + diff --git a/rksimpleVirus2.m b/rksimpleVirus2.m new file mode 100644 index 0000000..ca137ca --- /dev/null +++ b/rksimpleVirus2.m @@ -0,0 +1,47 @@ + +function rksimpleVirus2 + +r = 2.4; +a = 2; +b = 0.1; +c = 1; + + +x0 = r/a; + + + +xstar = r/a +vstar = b*r/(a*c) + +J = [-b c; -a*vstar r-a*xstar] +eig(J) +Delta = det(J) +tau = trace(J) + + + +y0 = [0.3509 1]; +tspan = [0 200]; + +[t,y] = ode45(@f5,tspan,y0); + +figure(1) +plot(t,y(:,1),t,y(:,2)) + +figure(2) +plot(y(:,1),y(:,2)) + + function yd = f5(t,y) + + yp(1) = c*y(2) - b*y(1); + yp(2) = y(2).*(r - a*y(1)); + + + + yd = [yp(1);yp(2)]; + + end % end f5 + +end + diff --git a/rossler.m b/rossler.m new file mode 100644 index 0000000..3ccd8b7 --- /dev/null +++ b/rossler.m @@ -0,0 +1,93 @@ +% function rossler + +function rossler + +a = 0.15; % 0.15 (defined phase) 0.3 (chaotic phase) +b = 0.4; % 0.4 +c = 8; % 6 8.0 + +y0 = [1 1 0]; +tspan = [1 30]; +[t,y] = ode45(@f5,tspan,y0); +[sz dum] = size(y); + +y0 =[y(sz,1) y(sz,2) y(sz,3)]; +tspan = [1 200]; +options = odeset('OutputFcn',@odephas3); +[t,y] = ode45(@f5,tspan,y0,options); + +figure(2) +plot(t,y(:,1),t,y(:,2),t,y(:,3)) + +figure(20) +colormap(jet) +plot3(y(:,1),y(:,2),y(:,3),'k') +set(gcf,'Color','white') +grid on + + +figure(3) +plot(y(:,1),y(:,3)) +title('x-z') + +figure(4) +plot(y(:,2),y(:,3)) +title('y-z') + +figure(5) +plot(y(:,1),y(:,2)) +title('x-y') + +PowS = 1; +% Power Spectrum +if PowS == 1 % Turn this off by setting PowS = 0 + [st dum] = size(t); + tshort = t(400:st,1); + yshort = y(400:st,1); + %yshort = u(400:st,1); + [st dum] = size(tshort); + dt = (max(tshort)-min(tshort))/st; + df = 1/st/dt; + fmax = 1/dt; + freq = df:df:fmax; + + yf = yshort(:,1); + F = fft(yf); + Pow = F.*conj(F); + + mid = round(st/4); + figure(6) + semilogx(freq(2:mid),Pow(2:mid)) + title('Power Spectrum') + xlabel('Frequency (Hz)') + + + [Pmax I] = max(Pow); + f0 = freq(I) + T0 = 1/f0 +% T0nat = 2*pi/sqrt(w02) +end + + + +d = y(:,1); +e = y(:,2); +f = y(:,3); +Printfile4('lout',t',d',e',f'); + +%keyboard + + + function yd = f5(t,y) + + yp(1) = -(y(2) + y(3)); + yp(2) = y(1) + a* y(2); + yp(3) = b + y(3)*(y(1) - c); + + yd = [yp(1);yp(2);yp(3)]; + + end % end f5 + + +end % end ltest + diff --git a/sinecircle1.m b/sinecircle1.m new file mode 100644 index 0000000..ab96eee --- /dev/null +++ b/sinecircle1.m @@ -0,0 +1,77 @@ +%sinecircle1.m + +clear + +%omega = 1.618033988749895; %exp(1)/2 pi/2 golden +%omega = 3/2; +omega = 0; +%omega = 1; +%omega = sqrt(2); +%omega = 21/13; +%omega = 13/8; +%omega = 8/5; +%omega = 5/3; +%omega = 2/3; + +testcase = 2; % 1 = iterate at a single g 2= scan g + +if testcase == 1 + gval = 80.; % 43 76 + glooplo = gval; + gloophi = gval; +elseif testcase == 2 % scan g + figure(3) + close + glooplo = 0; + gloophi = 2; +end + + +delg = 0.005; +for gloop = glooplo:delg:gloophi + + g =gloop; + + if testcase == 1 + disp(strcat('g = ',num2str(g))) + end + + thet0 = rand; + thet1last = thet0; + N = 100; + for loop = 1:N + + thet1emp = mod(thet1last + omega + g*sin(2*pi*thet1last),1); + thet1last = mod(thet1emp,1); + + end + + for loop = 1:N + + thet1(loop) = mod(thet1last + omega + g*sin(2*pi*thet1last),1); + thet1last = mod(thet1(loop),1); + + end + + if testcase == 1 + figure(1) + plot(thet1,thet2,'o') + axis([0 1 0 1]) + + figure(2) + plot(thet1,(thet2-thet1),'o') + axis([0 1 0 1]) + elseif testcase == 2 + x = g*ones(1,N); + y = thet1; + hold on + figure(3) + plot(x,y,'.k') + axis([0 gloophi 0 1]) + hold off + set(gcf,'Color','White') + end + +end + + diff --git a/zerodiag.m b/zerodiag.m new file mode 100644 index 0000000..765ac5b --- /dev/null +++ b/zerodiag.m @@ -0,0 +1,24 @@ +% [B, D] = zerodiag(A,altmin,varargin) +% Sets diagonal terms to zero +% D is the diagonal values + +function [B, D] = zerodiag(A,altmin,varargin) + +[sy,sx] = size(A); +B = A; +if sy~=sx + disp('in zerodiag not square'); + B = A; +else + if nargin == 1 + for loop = 1:sx + B(loop,loop) = 0; + D(loop) = A(loop,loop); + end + else + for loop = 1:sx + B(loop,loop) = altmin; + D(loop) = A(loop,loop); + end + end +end \ No newline at end of file