update 2/21/17

HW4/README.md

@@ -0,0 +1,37 @@
+# Homework #3
+## due 3/1/17 by 11:59pm
+
+
+1. Use your repository 'roots_and_optimization'. Document all the HW4 work under the
+heading `# Homework #4` in your `README.md` file
+
+    a. Create a function called 'collar_potential_energy' that computes the total
+    potential energy of a collar connected to a spring and sliding on a rod. As shown in
+    the figure given a position, xc, and angle, theta:
+
+    ![Collar-mass on an inclined rod](collar_mass.png)
+
+    The spring is unstretched when x_C=0.5 m. The potential energy due to gravity is:
+
+    PE_g=m x_C\*g\*sin(theta)
+
+    where m=0.5 kg, and g is the acceleration due to gravity, 
+
+    and the potential energy due to the spring is:
+
+    PE_s=1/2\*K \*(DL)^2
+
+    where DL = 0.5 - sqrt(0.5^2+(0.5-x_C)^2)
+
+    b. Use the `goldmin.m` function to solve for the minimum potential energy at xc when
+    theta=0. *create an anonymous function with `@(x) collar_potential_energy(x,theta)` in
+    the input for goldmin. Be sure to include the script that solves for xc*
+
+    c. Create a for-loop that solves for the minimum potential energy position, xc, at a
+    given angle, theta, for theta = 0..90 degrees. 
+
+    d. Include a plot of xc vs theta. `plot(theta,xc)` with
+
+    `![Steady-state position of collar on rod at angle theta](plot.png)`
+
+3. Commit your changes to your repository. Sync your local repository with github.

lecture_10/.ipynb_checkpoints/lecture_10-checkpoint.ipynb

@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

lecture_10/GaussNaive.m

@@ -0,0 +1,25 @@
+function [x,Aug] = GaussNaive(A,y)
+% GaussNaive: naive Gauss elimination
+%   x = GaussNaive(A,b): Gauss elimination without pivoting.
+% input:
+%   A = coefficient matrix
+%   y = right hand side vector
+% output:
+%   x = solution vector
+[m,n] = size(A);
+if m~=n, error('Matrix A must be square'); end
+nb = n+1;
+Aug = [A y];
+% forward elimination
+for k = 1:n-1
+  for i = k+1:n
+    factor = Aug(i,k)/Aug(k,k);
+    Aug(i,k:nb) = Aug(i,k:nb)-factor*Aug(k,k:nb);
+  end
+end
+% back substitution
+x = zeros(n,1);
+x(n) = Aug(n,nb)/Aug(n,n);
+for i = n-1:-1:1
+  x(i) = (Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);
+end

lecture_10/GaussPivot.m

@@ -0,0 +1,33 @@
+function [x,Aug,npivots] = GaussPivot(A,b)
+% GaussPivot: Gauss elimination pivoting
+%   x = GaussPivot(A,b): Gauss elimination with pivoting.
+% input:
+%   A = coefficient matrix
+%   b = right hand side vector
+% output:
+%   x = solution vector
+[m,n]=size(A);
+if m~=n, error('Matrix A must be square'); end
+nb=n+1;
+Aug=[A b];
+npivots=0; % initially no pivots used
+% forward elimination
+for k = 1:n-1
+  % partial pivoting
+  [big,i]=max(abs(Aug(k:n,k)));
+  ipr=i+k-1;
+  if ipr~=k
+    npivots=npivots+1; % if the max is not the current index ipr, pivot count
+    Aug([k,ipr],:)=Aug([ipr,k],:);
+  end
+  for i = k+1:n
+    factor=Aug(i,k)/Aug(k,k);
+    Aug(i,k:nb)=Aug(i,k:nb)-factor*Aug(k,k:nb);
+  end
+end
+% back substitution
+x=zeros(n,1);
+x(n)=Aug(n,nb)/Aug(n,n);
+for i = n-1:-1:1
+  x(i)=(Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);
+end

lecture_10/Tridiag.m

@@ -0,0 +1,22 @@
+function x = Tridiag(e,f,g,r)
+% Tridiag: Tridiagonal equation solver banded system
+%   x = Tridiag(e,f,g,r): Tridiagonal system solver.
+% input:
+%   e = subdiagonal vector
+%   f = diagonal vector
+%   g = superdiagonal vector
+%   r = right hand side vector
+% output:
+%   x = solution vector
+n=length(f);
+% forward elimination
+for k = 2:n
+  factor = e(k)/f(k-1);
+  f(k) = f(k) - factor*g(k-1);
+  r(k) = r(k) - factor*r(k-1);
+end
+% back substitution
+x(n) = r(n)/f(n);
+for k = n-1:-1:1
+  x(k) = (r(k)-g(k)*x(k+1))/f(k);
+end