%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIRST DAY MATLAB worksheet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 453-553
%%% this file helps to get familiar with basic elements of MATLAB
%%% for use in the Numerical Methods for PDE course
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% M. Peszynska,  2013/4/1
%% 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 MATLAB and just want to work on numerical PDEs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%					       part A 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% GETTING ACQUAINTED WITH MATLAB syntax, vectors etc. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
%% or you can type doc function
%%
f = inline('exp(2*t)');     
%%% you can also use
f = @(t)(exp(2*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('2*exp(2*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
%%%%%%%%%%%%%%%%%%%%
%%% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%                                         part B
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PLOTTING SOLUTIONS OF ELLIPTIC, HYPERBOLIC, AND PARABOLIC EQUATIONS: 
%%%%%%%%%%%%%% elliptic
x = linspace(-1,1); y = x;
[xx,yy]= meshgrid(x,y);
surfc(xx,yy,2+xx+yy); title('An affine function of x,y');
%% 
surfc(xx,yy,xx.^2+yy.^2); title('What is my Laplacian ?');
%% 
surfc(xx,yy,log(sqrt(xx.^2+yy.^2)); title('A fundamental solution of Laplace equation');
%%%%%%%%% some time dependent problems: at best, place these lines in an M-file
dt = .1; t = 0:dt:1;
%%% solution of first order hyperbolic equation
for n=1:length(t)
    plot(x,sin(x-t(n)));
    axis([-1 1 -1 1]);
    title(sprintf('Wave at t=%g',t(n)));
    pause;
end
%%% solution of second order hyperbolic equation u_tt=c^2 u_xx
x = linspace(0,1);
for n=1:length(t)
    plot(x,sin(pi*x)*cos(10*t);
    axis([-1 1 -1 1]);
    title(sprintf('Wave at t=%g   WHAT is c?',t(n)));
    pause;
end
%%% EXERCISE: come up with a plot of solution to u_tt=c^2(u_xx+u_yy)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% solution of a parabolic equation (diffusion=heat conduction) u_t=c^2 u_xx
x = linspace(0,1);
for n=1:length(t)
    plot(x,sin(pi*x)*exp(-pi*pi*t(n)));
    axis([0 1 -1 1]);
    title(sprintf('Temperature at t=%g... WHAT IS c?',t(n)));
    pause;
end
%%% EXERCISE: come up with a plot of solution to u_t=c^2(u_xx+u_yy)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% HOW TO SHOW ORDER OF ACCURACY/CONVERGENCE ORDER
%% Use functions error_loglog.m and error_table.m from Textbook website
%% to plot the error in approximating u''(x) by 1/h^2*(2u(x)-u(x-h)-u(x+h))
%% which is second order accurate.
%% Use u(x) = sin(pi*x) and check the value at x=1 when testing.
%%% Now derive and test a fourth order difference formula as in Exercise 1.2. 
%%%%%%%%%%%%%%%%%%%%
%% EXTRA credit: solve Exercises 1.2, 2.4, 2.5.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%