% Script 1 for Week 1 

% The very first code snippet below introduces the idea of a loop and fixed-point iteration, with the loop implemented by using the "up arrow and return." Be sure to point out the changing contents of variable x in the Workspace window as you iterate. Essentially all the programming ideas needed in the whole course are examined in the first week.

x=7;

% semicolon at the end of a line suppresses display of line output
% Look at the workspace;

x=sqrt(1+x)

% "New x" = "function of old x" so new x updates the value of previous x.
% Up arrow and return to repeat last command. Examine workspace.

% hit up arrow and return a bunch of times and see that the result is
% converging. To what? Anyone recognize this number?

format long
pi
format short
pi

% Watch in the workspace when you type the next few commands 

x=7,y=3
clear x
clear

x=pi:0.4:2*pi

H=sin(x)

figure(7)
plot(x,H)
figure(8)
plot(x,H,x,H,"r*")

% The following identifies the second, 4th-6th and last entry of H

H([2 4:6 8])

% The following squares those particular entries: notice "dot operation" which applies the operation to each of these specific entries, leaving the rest alone.

H([2 4:6 8])=H([2 4:6 8]).^2

% Notice error below, with "un-dotted" squaring operation. Always use dotted if you want the operation done on each entry.

H([2 4:6 8])^2

7./H([2 4:6 8])

% Now, your first "automatic" loop
% The pause command will cause execution to halt till you hit any key, rather than passing through the loop 4 times all at once.

x = 7
for k = 1:4
x = sqrt(1 + x)
pause
end


% This is called the "fixed point method" for solving an equation, 
% in this case we are finding a root of x-sqrt(1+x).

phi=x
x=-1:0.02:4

% Here is a graphical method for estimating phi, also called the "golden ratio."

y1 = x;
y2 = sqrt(1+x);
plot(x,y1,x,y2,phi,phi,'o')

% Logical operations or "|" , and "&"

b = 1;
a = 20;
(b ~= 0) & (a/b > 18.5)
(b ~= 0) & (a/b < 18.5)

% Difference between &, |, || and &&
% Run these groups of lines one at a  time.
% || and && are called the "short-circuit" "or" and "and" operators. The expression is evaluated as soon as it is determined with a short-circuit operator, which in this case is immediately after the first expression is determined to be false.

w=1:10^8;

tic
(b ~= 1) && (sum(w)==7)
toc

tic
(b ~= 1) & (sum(w)==7)
toc

% Below, both expressions are true. One is faster to evaluate.

tic
(b == 1) || (sum(w)==7)
toc


tic
(b == 1) | (sum(w)==7)
toc

% You might want to use the non-short circuit version of a logical expression if you wanted to make sure both expressions could actually be evaluated, possibly as an "error check" of some kind.

(1 == 1) || (1/[7 3]==7)

(1 == 1) | (1/[7 3]==7)

% Below, we use the while expression to calculate the same loop as above. We also implement a "counter" to keep track of how many times the loop is traversed. Do you think we have found the EXACT place where x=sqrt(1+x)? If not, why did the loop stop?

x=7
k=0
while x ~= sqrt(1+x)
k=k+1;
x = sqrt(1+x);
end
k
format long
x
format short

% Every number x has an "epsilon." Any number smaller than this epsilon, when added to x, is not distinguishable from x by MATLAB. This has to do with the way MATLAB stores numbers in memory. 

eps(7)
eps(14)
7==7+0.5*eps(7)

% MATAB thinks these two numbers are the same. But MATLAB knows the following two are different.

7==7+1.1*eps(7)

% Here is another loop that calculates phi to many places, with a step counter built-in.

x=7;
y = 0;
k=0;
while abs(x-y) > eps(x)
k=k+1;
y = x;
x = sqrt(1+x);
end
x
k

% Here is a very important "branching method" that causes MATLAB to do different things depending on a specific condition. Change the value of s in various ways and see what value MATLAB assigns to x. The first condition starts with if. The last starts with else. Intermediate branchings start with elseif. "elseif" and "else" are optional but "end" is not.

s=30;

if s>4
x=9
elseif s==4
x=10
else 
x=11
end

% In a loop increments need not be integers! Below we define two vectors and graph them as x,y pairs. In the context of the "for" command t takes on the 5 values 0, 0.02, 0.04, 0.06, 0.08 one after the other which are used in the calculation of the vector components x(i) and y(i). At each step i is incremented by 1. At the start of 6th step we add a 5th 0.2 and get 1.0, which is beyond 0.08, so we immediately exit (end) the loop without assigning that last value to t or doing anything else.

clear
i = 0;
for t = 0:0.02:0.08
i = i + 1;
x(i)=t;
y(i) = cos(t);
end

plot(x,y)

% Of course, as we saw above the same graph could be produced more easily using built-in vector methods. We would need to use the code above if subsequent values of x(i) or y(i) required knowledge of earlier values.

% Look at the workspace before and after typing the following commands. What changes?

t = 0:0.02:0.08
z=cos(t)

plot(t,z)

% For the instructor: you should go over a quick intro to the IDEA of a differential equation y'=f(t,y,etc) and Euler's method. Do not get bogged down, don't calculate ANY non-obvious analytical solutions, just the idea of a DE!! (You should mention that many---but by no means all---DEs have analytical solutions, and DE classes spend most of their time on methods to find these.) Maybe show a picture of a slope field and drawn-in solution curves, the graphical equivalent of Euler's method. 

% http://susanka.org/Classes/240INFO/01.PDF
% http://susanka.org/Classes/240INFO/A.PDF
% http://susanka.org/Classes/240INFO/C.PDF

% We will use MATLAB to do any messy arithmetic at the end of today's script.

% Here is a mathematical modeling example, for moderate velocities, for velocity of a falling body in air near Earth surface

% The differential equation for a falling object with velocity-dependent drag is
 
% vprime = g - cd/m*v^2 

% where cd is the drag coefficient and m is the mass and g is the gravitational acceleration. 

% This is a typical example that Chapra returns to often. DO NOT delve into the theory of this DE, or the physics associated with it beyond understanding what the terms mean. Chapra contains hundreds of examples like this and whole classes could be devoted to each one. We expect that the experts know their business and we get on with ours: the numerical implications.

% Carefully follow the workspace changes as you type in the following.


m=68.1
cd=0.25
g=9.81

t=12;

% We happen to know the exact solution at time as produced in a DE class, so we can compare our approximation technique to that and see how it behaves. At time t=12 we have:

vsol=sqrt(g*m/cd)*tanh(sqrt(g*cd/m)*t)

% We are going to approximate one "step" of a solution using Euler's method, as described in the text, starting with velocity zero. 

deltat=0.2
vapprox(1)=0
vapprox(2)=vapprox(1)+(g-cd/m*vapprox(1)^2)*deltat

% vapprox(2) is the (crude) approximation at time t=0.2.

% Now we'll look at 3 approximations which should be getting better and better because the time increment is getting smaller and smaller---each over an interval of twelve seconds.


deltat1=1;
vapprox1(1)=0;
for k=2:(12/deltat1+1);
vapprox1(k)=vapprox1(k-1)+(g-cd/m*vapprox1(k-1)^2)*deltat1;
end
t1=0:deltat1:12;


deltat2=0.5;
vapprox2(1)=0;
for k=2:(12/deltat2+1);
vapprox2(k)=vapprox2(k-1)+(g-cd/m*vapprox2(k-1)^2)*deltat2;
end
t2=0:deltat2:12;


deltat3=0.2;
vapprox3(1)=0;
for k=2:(12/deltat3+1);
vapprox3(k)=vapprox3(k-1)+(g-cd/m*vapprox3(k-1)^2)*deltat3;
end
t3=0:deltat3:12;

% Here is the exact solution for comparison. (Normally we won't have this, but since we do it is interesting to compare exact to approximation.)

vsol = sqrt(g*m/cd)*tanh(sqrt(g*cd/m)*t3);

plot(t3,vsol,'b-',t1,vapprox1,'o',t2,vapprox2,'*',t3,vapprox3,'d')

clf

plot(t3,vsol,t1,vapprox1,t2,vapprox2,t3,vapprox3)


% How close is close enough? And if you didn't have the exact solution how would you gauge how far away you were from the real solution?

% And if the velocities are close, what does that mean about the error in position estimate? As the sky-diver, wouldn't that be the real issue for you?

p1(1)=0;
for k=2:(12/deltat1+1);
p1(k)=p1(k-1)+vapprox1(k-1)*deltat1;
end

p2(1)=0;
for k=2:(12/deltat2+1);
p2(k)=p2(k-1)+vapprox2(k-1)*deltat2;
end

p3(1)=0;
for k=2:(12/deltat3+1);
p3(k)=p3(k-1)+vapprox3(k-1)*deltat3;
end

plot(t1,p1,'-',t2,p2,'-',t3,p3,'-')

% Would YOU rely on this estimate if you were a sky-diver and counting on its accuracy for YOUR life????