Review the Predator Prey notes. This assignment is to solve the Lotka-Volterra equations using the *backward* Euler method. Follow the procedure I used in class for the FHN equations. You may use the backward EulerND and newtonND examples from class. You will need to write a function called IvBE that will be passed through backward EulerND and into newtonND to evaluate the function value vector and Jacobian matrix needed by newtonND. It will be analogous to fhnBE from class. Also write a script that sets parameter values, timestep increment, etc., then calls backward EulerND and plots the solution. It would be wise to first solve the system using rk4 or eulerND to be sure you are getting a correct solution from backward EulerND. Call your script runLV. Use a time step of 0.001. Set the initial condition for chipmunk population, c, to 1.0 and the initial kitty population, k, to 0.1. At first, use the parameters: alpha = 1 beta = 10 delta=1 gamma=1 This is actually a fixed point of the system of ODEs, so the two populations will stay perfectly balanced at the initial values. Next, make the cats less hungry by reducing beta, then more hungry by increasing beta. Observe what happens to the two populations in the two cases. Then set beta back to 10 and make the cats longer lived by increasing gamma to 10. Again, observe what happens to the populations. Feel free to play with other parameter combinations. Upload runLV.m, LVBE.m and a plot of your solution for the parameters alpha=1, beta=15, delta=1, gamma=1. Your code must be appropriately commented and indented. function [F, J] = fhnBE (Unp1, paramsNewton) % function [F, J] = fhnBE (Unp1, params Newton) % % Used to solve the FHN system of ODEs using the multidimensional backward % Euler method. It is called by newtonND which is called by % backward EulerND. The system is % % % du/dt = c1*u* (u-a)*(1-u) - c2*u*v + Ileak. dv/dt b* (u-d*v) % % Applying the backward Euler method produces two nonlinear expressions % that must be solved in each timestep. % % let unp1 be u(n+1), the value of u in the next timestep we are trying to % find. similarly, vnp1 is v(n+1). un and vn are the known values of u and % v in the current timestep, u(n) and v(n). The two expressions are: % % f(unp1, vnp1) = un + (c1* unp1* (unp1-a)*(1-unp1) - c2unp1*vnp1 + Il )*dt unp1 % g(unp1, vnp1) = vn + b* (unp1 - d*vnp1) *dt % % The input vector Unp1 contains unp1 and vnp1. paramsNewton contains the % constants in the ODE as well as un, vn, and dt. It was packed up by % backward Euler ND and passed into newtonND, which passed it here. The % function returns F, which is f and g stacked into a column vector, and J, % the Jacobian matrix of F % unpack current estimates for unp1 and vnp1. The very first time % newton ND calls this function, these will be un and vn. newton ND will % improve these estimates iteratively until it finds values that when % plugged into f and g make both expressions 0. unp1 = Unp1 (1); vnp1 = Unp1(2); % unpack the constants in the ODE. These were passed into % backward Euler ND as params. a = params Newton (1); b = paramsNewton (2); c1= paramsNewton (3); c2 paramsNewton (4); d = paramsNewton (5); Il = paramsNewton (6); % unpack un and vn un vn params Newton (7); paramsNewton (8); % unpack t(n+1). Time does not appear in the RHSs of the FHN system, so % we don't need it. tnp1 = params Newton (9); %unused % unpack dt dt = paramsNewton (10); % we now have all the info we need to evaluate f and g. f = un + (c1*unp1* (unp1-a)*(1-unp1) - c2*unp1*vnp1 +Il) *dt g = vn + b* (unp1-d*vnp1) *dt - vnp1; F = [f;g]; % stack into column vector end % compute the elements of the Jacobian matrix. dfdunp1 is the % derivative of f with respect to unp1 - dfdunp1 c1*dt* ( -unp1* (unp1-a) + unp1*(1-unp1) + (unp1-a)*(1-unp1) ) - c2*dt*vnp1-1; % the other components are simpler... dfdvnp1 = -c2*dt*unp1; dgdunp1 b*dt; dgdvnp1 = -b*d*dt - 1; = % deriv of f wrt vnp1 % deriv of g wrt unp1 % deriv of g wrt vnp1 J = [dfdunp1 dfdvnp1; dgdunp1 dgdvnp1]; unp1; % assemble Jacobian. derivs of f in top row, of g in second row. Look % at the nonlinear 2D notes to see see why F and J are structured this % way. function [U, t] = backward Euler ND (UO, dt, tend, params, fcn) % function [U, t] = backward Euler ND (U0, dt, tend, params, fcn) % % Solves a *system* of ODEs using the backward Euler method. The calling % sequence is simlar to all our other ODE integrating functions. This uses % the same idea as backward Euler for a *single* ODE, but taking a backward % Euler timestep for a *system* of ODEs results in a *system* of nonlinear % algebraic expressions whose root must be found in each timestep. We use % the multidimensional Newton method to do this. For the FHN example, fcn % will be a handle to fhnBE, which evaluates these expressions % initialize solution array U and time vector t as before. U(1, :) = U0; t = 0:dt: tend; U(length(t), :) = 0; % timestep loop for n-1: length(t)-1 end end % Our initial guess for the solution U at timestep n+1 is just the % solution at timestep n. NewtonND expects its intitial guess for % the root to be in a column vector, so we transpose. guess = U(n, :)'; tol=1.e-8; % convergence tolerance for Newton ND % The params vector passed into this function contains the % constants in the ODEs. For the FHN system, this is a, b, c1, c2, d, % and Ileak. fcn, the function that evaluates the system of % nonlinear expressions we are finding the root to (and its % Jacobian matrix) also needs some other information: the current % U, sometimes t(n+1) and the timestep. We pack all this stuff up % into another parameter vector that we pass into newtonND. paramsNewton = [params U(n,:) t(n+1) dt]; % call newtonND. It will return a column vector that contains the % solution for the n+1 timestep. Transpose it and store it in U. U(n+1, ) = newton ND (guess, tol, fcn, paramsNewton)'; function [root, n] = newton ND (guess, tol, fcn, params) % function [Root, n] = newtonND (guess, tol, fcn, params) % % Find a root of a system of N nonlinear expressions using the % multidimensiuonal Newton method. Works for any N. Guess is % a column vector (length N) with an inital estimate for the root. % tol is convergence % tolerance. params contains constants to be passed to fcn. fcn is handle % to a MATLAB function that computes the nonlinear expressions and the % associated Jacobian. Output root is column vector with N elements. n is % the iteration count. % % fcn has form [F, J] = fcn (X, params). X is column vector with N elements % containing independent variables. params contains constants. F is a % column vector with N elements containing the expressions evaluted with X. % J is the NxN Jacobian matrix. % copy the guess to X. As the algorithm runs, X will contain our % current best estimate to a root. We should really be checking that % guess is a column vector before proceeding X = guess; % call fcn to evaluate the N expressions and the NXN Jacobian matrix. [F, J] = fcn (X, params); n=1; % initialize iteration counter % iterate until the F is within tol of 0. This is a little more % complicated than the 1D newton method because F is a vector that % contains the values of N different functions. The usual way to handle % this is to use the Euclidean length of the vector, which we can get % from the norm() function. For example, for a 2D system with F = [f;g], % norm (F) = sqrt (f^2 + g^2). while norm (F) > tol % the following line is meat of method. Refer to notes. We solve a % linear system to get corrections to X. dx will be a column vector % with N elements dX = J\(-F); % update X by adding the corrections we just calculated X = X+dX; % reevaluate the fucntions and Jacobian [F, J] = fcn (X, params); n=n+1; % increment counter end % we get here when X has converged to a root. Copy it to root and % return root = X;