%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIRST DAY MATLAB worksheet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 452-552
%%% this file helps to get familiar with basic elements of MATLAB
%%% for use in the Numerical Methods for ODE course
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% M. Peszynska,  2008/1/14
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% Please place yourself in one of the two categories
%%% depending on your knowledge of MATLAB:
%%% A) I have never used MATLAB or  I need  a review
%%% B) I am familiar with (at least) some of MATLAB syntax and commands
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% INSTRUCTIONS: please type every line yourself and observe its output.
%%%               Ignore all lines and text beginning with %
%%%               Try to solve exercises or leave them for later.
%%% If you selected A, start from the begining and then go through B. 
%%% Otherwise skip to B or review part A to see how much you remember.
%%% At the end you will find supplementary exercises to explore after class.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%					       part A 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1 + 3
a = 1
b = 3
5 * a + b^2 + 1/a
1
1;
1,
%%%% vectors and matrices
x = [1 2 3]
x = [1;2;3]
%%%
A = [0 1; -1 0]
%% 
%% vector and matrix norms
%%
norm(x,inf)
norm(x,1)
%
norm(A)
norm(A,inf)
norm(A,2)
%%% seeking eigenvalues and eigenvalue decomposition
eig(A)
help eig
[v,d] = eig(A)
v*d*inv(v)
v*d*v'
%%
%%%%%%%%%%%%%%%%%%%%%% 
%%% how to set up and plot functions (you can also use M-files)
%% 
help function
%%
f = inline('exp(t)');
h = 1e-1;
x = 1/pi;
f(x+h)-f(x)
(f(x+h)-f(x))/h
d = (f(x+h)-f(x))/h
%% compare true derivative with the approximation
fprime = inline('exp(t)');
d - fprime(x),
%% compare it for many values of h: set up h as a vector
h = [1e-1 1e-2 1e-3]
for j=1:size(h,2) d = (f(x+h)-f(x))/h,aer = abs(d - fprime(x)),pause,end
%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% EXERCISE:
%%% set up to be the function required in HW
%%% you can collect errors in a vector by modifying aer to be aer(j) 
%% then you can plot it against h (use plot(h,aer))
%%%%%%%%%%%%%% EXERCISE: 
%%% write an M-file square.m that computes ^2: call it logistic(x)
%%% here is another way to set up a function: 
logistic = inline('u.*(4-u)');
%%% plot this function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% learn how to plot
u = 0:.1:1 %% note:	 	%this makes y into a vector
plot(u,logistic(u))
%%% plot it along u^2 (quadratic monomial)
%%% try this:
u^2
%%% better:
u.^2
%%% plot the two together
plot(u,u.^2,'b',u,logistic(u),'r');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%					       part B
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% start with EXERCISE: 
%%% plot the logistic function and u^2 in various intervals 
%%% so that the graph is smooth (how to achieve it ?)
%%% to help determine Lipschitz constant as in HW 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% differential equations
%%%
%%% learn how to solve a differential equation symbolically (analytically)
%%%
help dsolve
dsolve('Dy = -y + sin(t) ')
yfun = dsolve('Dy = -y + sin(t) ')
pretty(yfun)
yfun = dsolve('Dy = -y + sin(t) ','y(0)=1')
pretty(yfun)
%%%%
%%%% how to plot direction fields (by hand, you can also use program dsolve)
%%%%
help quiver
[t,y] = meshgrid(0:.1:2,-1:.1:1);
funinline = inline('-y + sin(t)','t','y'); 
u = ones(size(t)); v = funinline(t,y);
quiver(t,y,u,v)
%%%
%%% numerical solution of ODEs  
%%% there are many functions: here we use ode45
%%%
help ode45
%%% you have to code funfile as M-file here or use funinline without the @ character
%%% for example ode45(funinline,[0 2],[1]);
ode45(@funfile,[0,2],[1]);
ode45(@funfile,[0,2],[0]);	%% changed initial conditions
ode45(@funfile,[0,10],[1]);     %% solve and plot in larger interval
%%%
%%% how to deal with systems: solve the system form HW numerically: 
%%% set it up as system y'=A*y, A shown below
%%% 
%%%
harmonic = inline('[0 1; -1 0] *y','t','y');
[t,y]=ode45(harmonic,[0,2*pi],[1 0]);
plot(t,y(:,1))
plot(t,y(:,2))
plot(y(:,1),y(:,2))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% additional material
%%%%%%%%%%%%%%%%%%%%
%%% representation of numbers and what is machine epsilon
%%%%%%%%%%%%%%%%%%%%
format long		%% to see more digits
pi
-2*pi+100+pi-100 	%% test if associative law applies
format short
-2*pi+100+pi-100 
eps			%% this is the machine accuracy
1+eps			%% that is, the smallest positive number added to 1 
			%% which makes the difference: 
			%% but you can see only that difference in ...
format hex
1 + eps
1 
1 + eps/2               %% eps/2 does not change the value
format
%%%%%%%%%%%%%% EXERCISE: 
%%% find the smallest representable number: find n such that 1/n| is not zero
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%