%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Example Octave/MATLAB programs to produce the graphs %% for Math 3520, HW 4, Spring 2026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Ex.2(c) % Plot the flow field in phase space for a gradient system % X'= - grad V(X) % with potential function V, including the level curves of V % Initialize the grid in phase space where arrows will % indicate the flow X'=(x',y') at position X=(x,y): px=-2:0.2:2; % x grid range and spacing py=-1:0.2:1; % y grid range and spacing [x,y]=meshgrid(px,py); % 2-D grid ral = 0.4; % Relative arrow length % Define the flow field and potential function dx = y+2*x.*y; dy = x+x.*x-y.*y; dd = sqrt(dx.*dx + dy.*dy) + 0.01; % Raw speeds, nonzero dy = dy./dd; dx = dx./dd; % Normalized unit arrows V= y.^3 -x.*y -x.*x.*y; % potential function, so X'= -grad V(X) clf; % clear any previous graphic quiver(x,y,dx,dy,ral); % Phase portrait = flow field hold on; % superimpose the next graphic contour(x,y,V,20); % contour plot = level curves title("Exercise 2(c)"); %%%% Ex.4 % Plot the flow field in phase space in order to identify % a trapping region for use with the Poincare-Bendixson theorem % Initialize the grid in phase space where arrows will % indicate the flow X'=(x',y') at position X=(x,y): px=-3:0.5:3; % x grid range and spacing py=-3:0.5:3; % y grid range and spacing [x,y]=meshgrid(px,py); % 2-D grid ral = 0.4; % Relative arrow length % Define the flow field and potential function dx = x-y-x.^3; dy = x+y-y.^3; dd = sqrt(dx.*dx + dy.*dy) + 0.01; % Raw speeds, nonzero dy = dy./dd; dx = dx./dd; % Normalized unit arrows clf; % clear any previous graphic quiver(x,y,dx,dy,ral); % Phase portrait = flow field title("Exercise 4"); % Superimpose the boundaries of trapping region R hold on pRo = -2:0.1:2; % outer boundary parameterization Outer = 2*ones(size(pRo)); plot(pRo, Outer); % Top plot(pRo, -Outer); % Bottom plot(Outer, pRo); % Right plot(-Outer, pRo); % Left pRi = 0:0.1:1/2; % inner boundary parameterization Inner = (1/2)*ones(size(pRi)); plot(pRi, Inner-pRi); % NE plot(pRi, pRi-Inner); % SE plot(-pRi, Inner-pRi); % NW plot(-pRi, pRi-Inner); % SW %%%% Ex. 6 a=-0.9:0.05:.9; pmu=2./(1-a.*a); nmu=-pmu; plot(a,nmu); hold on; plot(a,pmu); xlabel("a"); ylabel("mu"); title("Exercise 6: (a,mu) regions with Hopf bifurcation boundaries."); %%%% Ex.7 % Initialize the grid in phase space where arrows will % indicate the flow X'=(x',y') at position X=(x,y): px=-3:0.5:3; % x grid range and spacing py=-3:0.5:3; % y grid range and spacing [x,y]=meshgrid(px,py); % 2-D grid ral = 0.4; % Relative arrow length % mu < 0 % Define the flow field and potential function mu = -2; dx = mu*x+y-x.^3; dy = mu*y-x-2*y.^3; dd = sqrt(dx.*dx + dy.*dy) + 0.01; % Raw speeds, nonzero dy = dy./dd; dx = dx./dd; % Normalized unit arrows quiver(x,y,dx,dy,ral); % Phase portrait = flow field title("Exercise 7, mu = -2"); % mu = 0 % Define the flow field and potential function dx = y-x.^3; dy = -x-2*y.^3; dd = sqrt(dx.*dx + dy.*dy) + 0.01; % Raw speeds, nonzero dy = dy./dd; dx = dx./dd; % Normalized unit arrows quiver(x,y,dx,dy,ral); % Phase portrait = flow field title("Exercise 7, mu = 0"); % mu > 0 % Define the flow field and potential function mu = 2; dx = mu*x+y-x.^3; dy = mu*y-x-2*y.^3; dd = sqrt(dx.*dx + dy.*dy) + 0.01; % Raw speeds, nonzero dy = dy./dd; dx = dx./dd; % Normalized unit arrows quiver(x,y,dx,dy,ral); % Phase portrait = flow field title("Exercise 7, mu = +2"); %%%% Ex.8 %%%% Runge-Kutta 4th order solver: function X = rk4a(X0,F,steps,h) % Input % X0 = [x0,y0] Row vector, initial position % F = function handle, computes Xdot(x,y) % steps = number of time steps to compute (e.g., 100) % h = time increment per step (e.g., 1/steps ) % Output % X = (1+steps)x(2) matrix of (x,y) values along the % trajectory starting at X0 X = zeros(1+steps,2); % Solution to be computed X(1,:) = X0; % Initial point for n=1:steps f1 = F(X(n,:)); f2 = F(X(n,:)+f1*h/2); f3 = F(X(n,:)+f2*h/2); f4 = F(X(n,:)+f3*h); X(n+1,:) = X(n,:) +h*(f1+2*f2+2*f3+f4)/6; end end %% Ex.8(b) Some trajectories from various (x0,y0) % Write the Xdot function for use in rk4a() mu=1.1; % slightly greater than mu_c=1 Xdot = @(V) [V(1).*(1-V(1).^2-V(2).^2)+V(2).*(mu-V(2)./sqrt(V(1).^2+V(2).^2)), V(2).*(1-V(1).^2-V(2).^2)-V(1).*(mu-V(2)./sqrt(V(1).^2+V(2).^2)) ]; h=0.01; % Use same time step for same accuracy hold on X = rk4a([ 0,-2], Xdot, 1000,h); plot(X(:,1),X(:,2)) X = rk4a([ 0, 2], Xdot, 1000,h); plot(X(:,1),X(:,2)) X = rk4a([-2, 0], Xdot, 1000,h); plot(X(:,1),X(:,2)) X = rk4a([ 2, 0], Xdot, 1000,h); plot(X(:,1),X(:,2)) X = rk4a([ 2, 2], Xdot, 1000,h); plot(X(:,1),X(:,2)) X = rk4a([-2,-2], Xdot, 1000,h); plot(X(:,1),X(:,2)) X = rk4a([ 2,-2], Xdot, 1000,h); plot(X(:,1),X(:,2)) X = rk4a([-2, 2], Xdot, 1000,h); plot(X(:,1),X(:,2)) X = rk4a([ 0, 0], Xdot, 1000,h); plot(X(:,1),X(:,2)); X = rk4a([-0.5, 0], Xdot, 1000,h); plot(X(:,1),X(:,2)) X = rk4a([ 0.5, 0], Xdot, 1000,h); plot(X(:,1),X(:,2)) X = rk4a([ 0, -0.5], Xdot, 1000,h); plot(X(:,1),X(:,2)) X = rk4a([ 0, 0.5], Xdot, 1000,h); plot(X(:,1),X(:,2)) title("Exercise 8(b), mu = 1.1 > mu_c and various (x0,y0)"); %% Ex.8(c) Estimate periods in terms of mu-mu_c x0=0; y0=1; % Always start at (x0,y0)=(0,1) h=0.01; % Use same time step for same accuracy % Try various values of mu>mu_c=1, with mu-mu_c small % (Begin with larger values, so shorter periods T(mu).) mu=1.5; Xdot = @(V) [V(1).*(1-V(1).^2-V(2).^2)+V(2).*(mu-V(2)./sqrt(V(1).^2+V(2).^2)), V(2).*(1-V(1).^2-V(2).^2)-V(1).*(mu-V(2)./sqrt(V(1).^2+V(2).^2)) ]; clf; X = rk4a([x0,y0], Xdot, 550,h); plot(X(:,1),X(:,2)); hold on; plot(zeros(5),-2:2); clf; X = rk4a([x0,y0], Xdot, 560,h); plot(X(:,1),X(:,2)); hold on; plot(zeros(5),-2:2); % 5.5 1 clf; plot(r, mu-sin(r)); hold on; plot(r,zeros(size(r))); xlabel("Cycle radius r"); ylabel("Factor (mu - sin(r))"); title("Exercise 9, mu > 1 [mu = 1.5]"); mu=1; % mu = 1 clf; plot(r, mu-sin(r)); hold on; plot(r,zeros(size(r))); xlabel("Cycle radius r"); ylabel("Factor (mu - sin(r))"); title("Exercise 9, mu = 1"); mu=0.5; % 0 < mu < 1 clf; plot(r, mu-sin(r)); hold on; plot(r,zeros(size(r))); xlabel("Cycle radius r"); ylabel("Factor (mu - sin(r))"); title("Exercise 9, 0 < mu < 1 [mu = 0.5]"); mu=0; % mu = 0 clf; plot(r, mu-sin(r)); hold on; plot(r,zeros(size(r))); xlabel("Cycle radius r"); ylabel("Factor (mu - sin(r))"); title("Exercise 9, mu = 0"); mu= -0.5; % -1 < mu < 0 clf; plot(r, mu-sin(r)); hold on; plot(r,zeros(size(r))); xlabel("Cycle radius r"); ylabel("Factor (mu - sin(r))"); title("Exercise 9, -1 < mu < 0 [mu = -0.5]"); mu= -1; % mu = -1 clf; plot(r, mu-sin(r)); hold on; plot(r,zeros(size(r))); xlabel("Cycle radius r"); ylabel("Factor (mu - sin(r))"); title("Exercise 9, mu = -1"); mu= -1.5; % mu < -1 clf; plot(r, mu-sin(r)); hold on; plot(r,zeros(size(r))); xlabel("Cycle radius r"); ylabel("Factor (mu - sin(r))"); title("Exercise 9, mu < -1 [mu = -1.5]"); %%%% Ex.10 % Plot the flow field in phase space in order to identify % a homoclinic bifurcation at mu_c ~ 0.66 % Initialize the grid in phase space where arrows will % indicate the flow X'=(x',y') at position X=(x,y): px=-3:0.5:3; % x grid range and spacing py=-3:0.5:3; % y grid range and spacing [x,y]=meshgrid(px,py); % 2-D grid ral = 0.4; % Relative arrow length muc = 0.066; % Critical mu. approximately % Define the flow field and potential function for mu < mu_c mu =-1; dx = mu*x+y-x.^2; dy = -x+mu*y+2*x.^2; dd = sqrt(dx.*dx + dy.*dy) + 0.01; % Raw speeds, nonzero dy = dy./dd; dx = dx./dd; % Normalized unit arrows clf; % clear any previous graphic quiver(x,y,dx,dy,ral); % Phase portrait = flow field title("Exercise 10, mu = -1 < mu_c"); % Define the flow field and potential function for mu = mu_c mu =0.066; dx = mu*x+y-x.^2; dy = -x+mu*y+2*x.^2; dd = sqrt(dx.*dx + dy.*dy) + 0.01; % Raw speeds, nonzero dy = dy./dd; dx = dx./dd; % Normalized unit arrows clf; % clear any previous graphic quiver(x,y,dx,dy,ral); % Phase portrait = flow field title("Exercise 10, mu = mu_c ~ 0.066"); % Define the flow field and potential function for mu > mu_c mu =2; dx = mu*x+y-x.^2; dy = -x+mu*y+2*x.^2; dd = sqrt(dx.*dx + dy.*dy) + 0.01; % Raw speeds, nonzero dy = dy./dd; dx = dx./dd; % Normalized unit arrows clf; % clear any previous graphic quiver(x,y,dx,dy,ral); % Phase portrait = flow field title("Exercise 10, mu = 2 > mu_c");