1. Describe the system of differential equations to be solved numerically. For this, save the following function as a file gain.m : ------------------------- CUT HERE ----------------------------- function xdot=gain(t,x) % This function codes the equations for the % model of nonlinear gain control in the % retina (Example 6, p. 37) tau_a=1; tau_b=1; L=1; B=x(1); A=x(2); Bdot=(-B + L/(1+A))/tau_b; Adot=(-A+2*B)/tau_a; xdot=[Bdot; Adot]; -----------------------------CUT HERE ---------------------------- 2. Run Matlab: matlab in the MATLAB command window type the following commands (skip the comments following %) and press ENTER at the end of each line. diary save_your_work % saves your work to file 'save_your_work' in % current directory clear all % clear the work space tspan=[0 50]; % defines the time interval for integration x0=[1; 0.5]; % defines initial conditions [t, x]=ode23('gain', tspan, x0); % calls MATLAB function ode23 for % numerical integration of the ODE % (described in gain.m) with initial data in x0 on the % time interval defined in tspan % The results are stored in t and x % t is a vector of discrete times % for each value of t, x contains a pair of values (x1(t), x2(t)) % for position and velocity, respectively B=x(:,1); % extract the first column with the values of B A=x(:,2); % extract the second column for A figure(1) % open a new figure plot(t,B,'-b') % plot B versus t figure(2) % open another figure plot(t,A) % plot A versus t figure(3) % alternatively we can combine two plots plot(t,B,'-b') % plot B versus t hold on % plot(t,A,'--r') % plot x2 vs t xlabel('time, t') % supply lablels ylabel('B (-), and A (--)') print -deps batraces % save this figure to file batraces.eps diary off % alternatively one can use the menu on the top of the matlab figure % to save or print the figure. % For this, choose 'File' and 'Save as', and choose from the list % of available formats. % 3. Save the following fragement as pln.m This code plots a phase portrait for the present problem. In matlab, type 'pln' and press 'enter'. ------------------------------BEGIN CUT ------------------------------------ figure(4) % open another figure hold on tspan=[0 20]; for i=0.2:0.2:1.8 x0=[0; i]; % defines initial conditions [t, x]=ode23('gain', tspan, x0); plot(x(:,1), x(:,2)) drawnow [t, x]=ode23('gain', tspan, [i; 0]); plot(x(:,1), x(:,2), '-b') [t, x]=ode23('gain', tspan, [i; 2]); plot(x(:,1), x(:,2), '-b') [t, x]=ode23('gain', tspan, [2; i]); plot(x(:,1), x(:,2), '-b') drawnow end axis([0.2 1.0 0.4 1.8]) % Next, plot the nullclines x1=0:0.05:2; y1=2*x1; % A-nullcline plot( x1, y1, '-r', 'linewidth', 2) % plots A-nullcline plot( 1./(1+x1), x1, '-g', 'linewidth', 2) % plots B-nullcline; note the dot % following in the expression of 1./(1+x1) % It indicates that the division in arrays % should be executed in term-by-term manner % Compute the coorinates of the fixed point in the first quadrant % b_bar=(-1+sqrt(1+8*1))/4; a_bar=2*b_bar; plot(b_bar, a_bar, 'ok', 'linewidth', 3) % indicate the fixed point on the plot % Compute the Jacobian: DF=[-1 -1/(1+a_bar)^2; 2 -1]; % ; indicates end of row %Find the eigenvalues and the eigenvectors of DF [V, D]=eig(DF); % eig is a matlab function which computes eigenvalues and eigenvectors % see http://www.mathworks.com/access/helpdesk/help/techdoc % /index.html?/access/helpdesk/help/techdoc/ref/eig.html % The eigenvalues are complex. Thus, the eigenvectors can be taken complex conjugate % We shall need the real and imaginary part of one of them v1=real(V(:,1)); v2=imag(V(:,1)); % extract the eigenvectors v1=v1/norm(v1); v2=v2/norm(v2); % normalize s=-1:1; plot(b_bar+v1(1)*s, a_bar+v1(2)*s, '--k', 'linewidth', 1) % In the present case the subspaces plot(b_bar+v2(1)*s, a_bar+v2(2)*s, '--k', 'linewidth', 1) % spanned by the eigenvectors are % are not as informative as they would be % if the eigenvalues were real. But we plot them % anyway. % Finally, add the axis lables and the title xlabel('B') ylabel('A') title('Phase Plane') print -deps pplane % save the figure -------------------------------------------END CUT--------- 4. Below, calculate the fixed point and the analyze the linearized system % Compute the coorinates of the fixed point in the first quadrant b_bar=(-1+sqrt(1+8*1)/4; a_bar=2*b_bar; % Compute the Jacobian: DF=[-1 -1/(1+a_bar)^2; 2 -1]; % ; indicates end of row % Find the eigenvalues and the eigenvectors of DF [V, D]=eig(DF) % eig is a matlab function which computes eigenvalues and eigenvectors % see http://www.mathworks.com/access/helpdesk/help/techdoc % /index.html?/access/helpdesk/help/techdoc/ref/eig.html % Extract the eigenvalues lambda_1=D(1,1) lambda_2=D(2,2) % and the eigenvectors v1=V(:,1) v2=V(:,2) v1=v1/norm(v1); v2=v2/norm(v2); % normalize