A6 Solutions

P1.

function up = SlopeF(t,u)
global nEvals

r = (u(1)^2 + u(3)^2)^(3/2);
up = zeros(4,1);
up = [u(2);-u(1)/r; u(4);-u(3)/r];
nEvals = nEvals + 1;



  function  [x,y,FunEvals] = ODE23Orbit(Tfinal,x0,xp0,y0,yp0,abserr)
% x and y are scalars that estimate (x(Tfinal),y(Tfinal)) where 
% x0 and xp0 are the x-coordinate and x-velocity at t = 0. 
% y0 and yp0 are the y-coordinate and y-velocity at t = 0.
% Uses ODE23 with the absolute error tolerance set to abserr.
% FunEvals is the number of required function evaluations.

u0 = [x0;xp0;y0;yp0];
global nEvals
nEvals = 0;
options = odeset('RelTol',abserr);
[t,ufinal] = ode23(@SlopeF,[0,Tfinal],u0,options);
x = ufinal(length(t),1);
y = ufinal(length(t),3);
FunEvals = nEvals;


  function  [x,y] = AB2Orbit(Tfinal,N,x0,xp0,y0,yp0)
% x and y are column N+1 vectors where (x(i),y(i)) is an
% estimate of the orbiting body's location at time t = (i-1)h
% where h = Tfinal/N. The second order Adams-Bashforth method
% is used with fixed stepsize h. x0 and xp0 are the x-coordinate
% and x-velocity at t = 0. y0 and yp0 are the y-coordinate and y-velocity
% at t = 0.

h = Tfinal/N;
x = zeros(N+1,1); x(1) = x0;
y = zeros(N+1,1); y(1) = y0;

t_c = 0;
u_c = [x0;xp0;y0;yp0];
f_c = SlopeF(t_c,u_c);

% First step by RK2

fr = SlopeF(t_c,u_c + h*f_c);
u_c = u_c + (h/2)*(f_c + fr);
x(2) = u_c(1); y(2) = u_c(2);
t_c = t_c + h;
f_minus = f_c;
f_c = SlopeF(t_c,u_c);

for i=3:N+1
    % Get x((i-1)h) and y((i-1)h))
    u_c = u_c + (h/2)*(3*f_c - f_minus);
    x(i) = u_c(1); y(i) = u_c(3);
    t_c = t_c + h;
    f_minus = f_c;
    f_c = SlopeF(t_c,u_c); 
end

% Grading scheme
%
%  -2 if not using Runge-Kutta for first step.

P2.

function u = LShape(N,F,G)

M = (N-1)*(2*N-1);         % Number of interior meshpoints in the top and bottom rectangles.
P = 2*M + (N-1);           % Total number of unknowns

% Set up the right hand side P-by-1 vector b

b = zeros(P,1);
h = 1/N; h_square = h*h;
y = 3;
k = 0;
for i=1:3*N-1
    % Process the ith horizontal line of meshpoints (starting from top)
    y = y-h % The y-value associated with the ith meshpoint line.
    x = 0;   % Reset x to the left edge value.
    
    if i<=2*N
        % Top Rectangle
        for j=1:N-1
            k=k+1;      
            x=x+h;      b(k) = -h_square*feval(F,x,y);
            if j==1 ,   b(k) = b(k) + feval(G,0,y); end    % West  neighbor on boundary
            if j==N-1,  b(k) = b(k) + feval(G,1,y); end    % East  neighbor on boundary
            if i==1,    b(k) = b(k) + feval(G,x,3); end    % North neighbor on boundary 
        end
        
    elseif i==2*N+1
        % The top meshpoint line in the bottom rectangle is special
        for j=1:2*N-1
            k=k+1; 
            x = x+h;      b(k) = -h_square*feval(F,x,y);
            if j==1,      b(k) = b(k) + feval(G,0,y); end % West  neighbor on the boundary
            if j==2*N-1,  b(k) = b(k) + feval(G,2,y); end % East  neighbor on the boundary
            if j>=N,      b(k) = b(k) + feval(G,x,1); end % North neighbor on the boundary
        end
        
    else
        % Bottom Rectangle
        for j=1:2*N-1
            k=k+1; 
            x = x+h;      b(k) = -h_square*feval(F,x,y);
            if j==1,      b(k) = b(k) + feval(G,0,y); end  % West  neighbor on boundary
            if j==2*N-1,  b(k) = b(k) + feval(G,2,y); end  % East  neighbor on boundary
            if i==3*N-1,  b(k) = b(k) + feval(G,x,0); end  % South neighbor on boundary
        end
    end
end

% Solve using FastPS. The key matrices...

%     A1 = L(2N-1,N-1)
%     A2 = eye(N-1,N-1) - diag(ones(N-2,1),1) - diag(ones(N-2),1),1)
%     A3 = L(N-1,2n-1)
%     C1 = E(:,M-N+2:M),  E = eye(M,M)
%     C2 = E(:,1:N-1),    E = eye(M,M)

% Compute A2_tilde = A2 - C1'*(A1^{-1}C1) - C2'*(A2^{-1}C2)
  
  A2_tilde = 4*eye(N-1,N-1) - diag(ones(N-2,1),1) - diag(ones(N-2,1),-1);
  
  for q=1:N-1
      w1 = zeros(M,1); w1(M-N+1+q) = 1; z1 = FastPS(w1,2*N-1,N-1); 
      w2 = zeros(M,1); w2(q) = 1;       z2 = FastPS(w2,N-1,2*N-1);
      A2_tilde(:,q) = A2_tilde(:,q) - z1(M-N+2:M) - z2(1:N-1);
  end

% Compute b2_tilde = b2 - C1'*(A1^{-1}b1) - C2'*(A3^{-1}b3) where

  b1 = b(1:M); b2 = b(M+1:M+N-1); b3 = b(M+N:P);

  z1 = FastPS(b1,2*N-1,N-1);
  z2 = FastPS(b3,N-1,2*N-1);
  b2_tilde = b2 + z1(M-N+2:M) + z2(1:N-1);

  % Solve A2_tilde *u2 = b2_tilde
  u2 = A2_tilde \ b2_tilde;

  b1_tilde = b1; b1_tilde(M-N+2:M) = b1_tilde(M-N+2:M) +  u2; u1 = FastPS(b1_tilde,2*N-1,N-1);
  b3_tilde = b3; b3_tilde(1:N-1)   = b3_tilde(1:N-1)   +  u2; u3 = FastPS(b3_tilde,N-1,2*N-1);
  
  u = [u1;u2;u3];
  
% Grading scheme
%
% -4 if space usage is more than linear in # of unknowns. Specifically,
%    only A2 needs to be formed explicitly. C1 and C2 should be generated one
%    column at a time as needed.
% -2 if setting up RHS incorrectly, e.g. incorrect boundary conditions.
%
%