PDE Qual

Old Qual from Purdue

Hermite Polynomials and Gaussian distribution Links

Gaussian functions from Wiki

Gaussian integrals from Wiki

Normal distribution from Wiki

Construct Hermite Polynomials from Gaussian distribution functions (see copy below)

Hermite polynomials from Wiki

Addition on Hermite polynomials


Construct Hermite Polynomials from Gaussian distribution functions, copied from link above:

Let us look at the Hermite polynomials in somewhat more detail. Consider a Gaussian function, .

.

We write

.

This defines the Hermite polynomials.

.

The parity of the Hermite polynomials is (-1)ⁿ.

H_n(z) is a nth degree polynomial in z.

Proof:

This statement holds for n = 1 and n = 2.

Let . Then

.

If H_n-1(z) is a polynomial of degree n-1, the H_n(z) is a polynomial of degree n. The statement therefore holds for all n.

The generating function

Consider F(z+l) = exp(-(z+l)²).

A Taylor series expansion, , yields

.

Let l’ = -l. Then

.

. (Taylor series expansion)

is called the generating function for the Hermite polynomials H_n(z).

Recurrence relations

We have already found

.

Other recurrence relations exist.

and combining the first two relations

.

[ is obtained by differentiating with respect to z,

, and equating terms of equal power of l.

is obtained by differentiating with respect to l.]

Summary

The Hermite polynomials H_n(z) are defined by

or .

The generating function of the Hermite polynomials is

.

The recurrence relations are

, ,

and

.

The differential equation satisfied by the Hermite polynomials is

.

The normalization and orthogonality condition are

Basic probability in Matlab

Mean value

A =
1 2 0
2 3 0

mean(A) = 1.5000 2.5000 0 % Mean of each column

mean(A, 2) = % Mean of each row
1.0000
1.6667


Variance

• V = var(X) returns

- for vectors X: its variance.
- for matrices: a row vector containing the variance of each column of X.

• V = var(X,w,dim)
  

takes the variance along the dimension dim of X.
Pass in 0 for w to use the default normalization by N-1, or 1 to use N.
The variance is the square of the standard deviation (STD)

Example:
A =
1 2 0
2 3 0

var(A) =
0.5000 0.5000 0 % variance of columns

var(A, 0, 2) = % variance of rows
1.0000
2.3333

Links for "fun"

Topics - Online mag. for English learners

List of traditional children games

Deep fun

Test strong convergence of Milstein Mthd for SDE

Below is the file

******************************************

%MILSTRONG  Test strong convergence of Milstein: vectorized
%
% Solves   dX = r*X*(K-X) dt + beta*X dW,  X(0) = Xzero,
%        where r = 2, K= 1, beta = 1 and Xzero = 0.5.
%
% Discretized Brownian path over [0,1] has dt = 2^(-11).
% Milstein uses timesteps 128*dt, 64*dt, 32*dt, 16*dt (also dt for reference).
%
% Examines strong convergence at T=1:  E | X_L - X(T) |.
% Code is vectorized: all paths computed simultaneously.

randn('state',100)
r = 2; K = 1; beta = 0.25; Xzero = 0.5;   % problem parameters
T = 1; N = 2^(11); dt = T/N;                       % from 0 to T=1. N = 2^11 steps, each of length dt            
M = 500;                                                            % number of paths sampled
R = [1; 16; 32; 64; 128];                               % Milstein stepsizes are R*dt, 5 R's for 5 paths

dW = sqrt(dt)*randn(M,N);                 % Brownian increments for each path each timestep
Xmil = zeros(M,5);                                  % preallocate array: M paths w. corresp. 5 R
for p = 1:5                                                   % Consider each time step separately
     Dt = R(p)*dt; L = N/R(p);                % Dt = new time step
                                                                       % L timesteps of size Dt = R dt
     Xtemp = Xzero*ones(M,1);            % All M paths begin at Xzero
     for j = 1:L                          % For each step
          Winc = sum(dW(:,R(p)*(j-1)+1:R(p)*j),2); 
          Xtemp = Xtemp + Dt*r*Xtemp.*(K-Xtemp) + beta*Xtemp.*Winc ...
              + 0.5*beta^2*Xtemp.*(Winc.^2 - Dt); % See below
end
     Xmil(:,p) = Xtemp;       % store Milstein solution at t = 1
end

Xref = Xmil(:,1);                                       % Reference solution: "exact" soln when timestep = dt 
Xerr = abs(Xmil(:,2:5) - repmat(Xref,1,4));   % Error in each path
mean(Xerr);                                              % Mean pathwise errors
Dtvals = dt*R(2:5);                                % Milstein timesteps used

subplot(224)                                            % lower RH picture
loglog(Dtvals,mean(Xerr),'b*-'), hold on
loglog(Dtvals,Dtvals,'r--'), hold off  % reference slope of 1
axis([1e-3 1e-1 1e-4 1])
xlabel('\Delta t')
ylabel('Sample average of | X(T) - X_L |')
title('milstrong.m','FontSize',10)

%%%% Least squares fit of error = C * Dt^q %%%%
A = [ones(4,1), log(Dtvals)]; rhs = log(mean(Xerr)');
sol = A\rhs; q = sol(2)
resid = norm(A*sol - rhs)

******************************************

Notes:

-

-

Colon (:) in Matlab

colon (:)

Create vectors, array subscripting, and for-loop iterators.

Description

The colon is one of the most useful operators in MATLAB. It can create vectors, subscript arrays, and specify for iterations.The colon operator uses the following rules to create regularly spaced vectors:

j:k is the same as [j,j+1,...,k] (empty if j > k)
j:i:k is the same as [j,j+i,j+2i, ...,k]
is empty if i == 0, if i > 0 and j> k, or if i <> where i, j, and k are all scalars.

Below are the definitions that govern the use of the colon to pick out selected rows, columns, and elements of vectors, matrices, and higher-dimensional arrays:

A(:,j) is the jth column of A
A(i,:) is the ith row of A
A(:,:) is the equivalent two-dimensional array. For matrices this is the same as A.
A(j:k) is A(j), A(j+1),...,A(k)
A(:,j:k) is A(:,j), A(:,j+1),...,A(:,k)
A(:,:,k) is the kth page of three-dimensional array A.
A(i,j,k,:) is a vector in four-dimensional array A. The vector includes A(i,j,k,1), A(i,j,k,2), A(i,j,k,3), and so on.
A(:) is all the elements of A, regarded as a single column. On the left side of an assignment statement, A(:) fills A, preserving its shape from before. In this case, the right side must contain
the same number of elements as A.

Examples

Using the colon with integers,
D = 1:4 results in D = 1 2 3 4

Using two colons to create a vector with arbitrary real increments between the elements,
E = 0:.1:.5 results in E = 0 0.1000 0.2000 0.3000 0.4000 0.5000

The command
A(:,:,2) = pascal(3)
generates a three-dimensional array whose first page is all zeros.

A(:,:,1) =
0 0 0
0 0 0
0 0 0

A(:,:,2) =
1 1 1
1 2 3
1 3 6

Brownian Path Simulation using Matlab

%BPATH1  Brownian path simulation

randn('state',100)           % set the state of randn to be 100
T = 1; N = 500; dt = T/N;    % Length = 1, do 500 steps, each step measures dt=1/500
dW = zeros(1,N);             % preallocate arrays ...
W = zeros(1,N);              % for efficiency

dW(1) = sqrt(dt)*randn;      % first approximation outside the loop ...
W(1) = dW(1);                % since W(0) = 0 is not allowed
for j = 2:N
 dW(j) = sqrt(dt)*randn;   % general increment, dW ~ sqrt(dt)
 W(j) = W(j-1) + dW(j);    % then go a step of dW(j) from W(j-1) to W(j)
end

plot([0:dt:T],[0,W],'r-')    % plot W against t: plot from time 0 to time T, increment dt
xlabel('t','FontSize',16)
ylabel('W(t)','FontSize',16,'Rotation',0)

Note: Matlab starts array from 1, not 0. So in order to include the starting point W(0) = 0 into the graph, we need to include it in the plotting command by doing plot([0:dt:T],[0,W])

Following is a more elegant and more efficient version of bpath1.m :
The "for" loop is replaced by a "vectorized" command
(KEY to writing efficient MATLAB)code (See DJHigham & NJHigha, MATLAB guide, SIAM, Philly 2000)


%BPATH2  Brownian path simulation: vectorized

randn('state',100)          % set the state of randn
T = 1; N = 500; dt = T/N;

dW = sqrt(dt)*randn(1,N);   % generate the array of increments
W = cumsum(dW);             % cumulative sum, W(j)= dW(1) + ... + dW(j)

plot([0:dt:T],[0,W],'r-')   % plot W against t, red line
xlabel('t','FontSize',16)
ylabel('W(t)','FontSize',16,'Rotation',0)

A higher simulation of Brownian path:
Average from about 1000 different paths. Sketch the average path and 5 of the brownian paths.
Here we evaluate the function u(W(t)) = exp(t + W(t)/2) along 1000 discretized Brownian paths.
The average of u(W(t)) over these paths is plotted with a solid blue line.
Five individual paths are also plotted using a dashed red line.

%BPATH3  Function along a Brownian path

randn('state',100)                                % set the state of randn
T = 1; N = 500; dt = T/N; t = [dt:dt:1];          % t = row vector, N components, from dt to 1, increments dt.

M = 1000;                                         % M paths simultaneously
dW = sqrt(dt)*randn(M,N);                         % increments for M paths. 1 path 1 row. Each row N increments ...
                                                % i.e., dW is an M-by-N array st dW(i,j) gives the increment dWj

W = cumsum(dW,2);                                 % cumulative sum accross the second dimension (column)
                                                % i.e., each row of W is the cumsum of the corresp. row of dW

U = exp(repmat(t,[M 1]) + 0.5*W);                 % repmat(t,[M 1]) form a "matrix" M row 1 column, each "entry" is t.
                                                % With t a N-dim row vector, this put M t's stack together,
                                                % forming a M by N matrix.
                                                % then add it to the matrix .5*W
                                                  % then take exp of this matrix
Umean = mean(U);                                  % compute columwise averages
plot([0,t],[1,Umean],'b-'), hold on               % plot mean over M paths
                                                 % x-axis is 0 and the vector t
                                                 % y-axis is 1. At time 0 the path begin at y=1. WHY???
plot([0,t],[ones(5,1),U(1:5,:)],'r--'), hold off  % plot 5 individual paths
                                                 % discrete "time" from 0 to 1
                                                 % over 5 paths, all begin at 1
                                                 % paths are first 5 rows of U
xlabel('t','FontSize',16)
ylabel('U(t)','FontSize',16,'Rotation',0,'HorizontalAlignment','right')
legend('mean of 1000 paths','5 individual paths',2)

averr = norm((Umean - exp(9*t/8)),'inf')          % sample error. Infinity norm of matrix
                                                 % largest row sum of A.
                                                 % w.o. 'inf' - 1-norm, largest column sum

To be explored....

Matlab - Basic Plotting

To plot one function in one graph (line plot)

t = 0 : pi/100 : 2*pi; %create a vector of values in the range [0, 2pi] in increments of pi/100
y = sin(t); %use this vector to evaluate the sine function over that range
plot(t,y)
grid on % Turn on grid lines for this plot

To plot several functions in one graph (several lines plot)

y = sin(t);
y2 = sin(t-0.25);
y3 = sin(t-0.5);
plot(t,y,t,y2,t,y3)

with lines type specified

t = 0:pi/100:2*pi;
y = sin(t);
y2 = sin(t-0.25);
y3 = sin(t-0.5);
plot(t,y,'-',t,y2,'--',t,y3,':')

Colors, Line Styles, and Markers

The basic plotting functions accepts character-string arguments that specify various line styles, marker symbols, and colors for each vector plotted. In the general form,

plot(x,y,'linestyle_marker_color')

linestyle_marker_color is a character string (delineated by single quotation marks) constructed from

1. A line style (e.g., dashed, dotted, etc.)
2. A marker type (e.g., x, *, o, etc.)
3. A predefined color specifier (c, m, y, k, r, g, b, w)

For example,

plot(x,y,':squarey')%plots a yellow dotted line and places square markers at each data point.

If you specify a marker type, but not a line style, MATLAB draws only the marker.
The specification can consist of one or none of each specifier in any order.

For example, the string'go--' defines a dashed line with circular markers, both colored green.

You can also specify the size of the marker and, for markers that are closed shapes, you can specify separately the colors of the edges and the face.

x = -pi:pi/10:pi; % x is from -pi to pi, with increments pi/10
y = tan(sin(x)) - sin(tan(x)); % definition of y
plot(x,y,'--rs',... % plot the graph y = f(x), a red dashed line with squared markers
'LineWidth',2,... % a line with width of two points
'MarkerEdgeColor','k',... % the edge of the marker colored BLACK
'MarkerFaceColor','g',... % the facr of the marker colored green
'MarkerSize',10) % the size of the marker set to 10 points

Note: the three dots '...' tells Matlab that something more will be added in the next line.
W.o. them, Matlab will interpret the commend to be ended right when the line ends.

Adding plot to an existing graph

semilogx(1:100,'+') % screate a semilogarithmic plot
hold all % hold plot and cycle line colors
plot(1:3:300,1:100,'--') % then add a linear plot
hold off % stop the holding process
grid on % Turn on grid lines for this plot

Plotting only the data points

For example, given two vectors

x = 0:pi/15:4*pi; % x from 1 to 4pi, with increments pi/15.
y = exp(2*cos(x)); % y is a function in x
plot(x,y, 'r+') % call plot with only a color and a marker specifier, here a red plus sigh at each data point.

Plotting Markers and Lines

x = 0:pi/15:4*pi;
y = exp(2*cos(x)); % x and y as above
plot(x,y,'-r',... % data as a red, solid line
x,y,'ok') % adds circular markers with black edges at each data point



Some options:
'a' - none
'b' = blue
'c' = cyan
'd' = looks like {}
'e' - none
'f' - none
'g'
'h' = star
'k' = black
'l' - none
'm' = pink/purple
'n' - none
'o' = circle
'p' = star
'q' - none
'r' = red
's' = square
't' - none
'u' - none
'v' - looks like gradient notation
'w' = white
'y' = yellow
'z' - none

'^' = upward triangle
'*' = snowflake
':' = tiny dashed
'<' = rightward triangle '>' = leftward triangle
'.' = tinny colored dot

Matlab - Arrays

zero

B = zeros(n) returns an n-by-n matrix of zeros. An error message appears if n is not a scalar.
B = zeros(m,n) or B= zeros([m n]) returns an m-by-n matrix of zeros.
B = zeros(m,n,p,...) or B= zeros([m n p ...]) returns an m-by-n-by-p-by-... arrayof zeros.
B = zeros(size(A)) returns an array the same size as A consisting of all zeros.

Remarks

The MATLAB language does not have a dimension statement;

MATLAB automatically allocates storage for matrices.

Nevertheless, for large matrices, MATLAB programs may execute faster if the zeros function is used to set aside storage for a matrix whose elements are to be generated one at a time,
or a row or column at a time.

For example

x = zeros(1,n);
for i = 1:n, x(i) = i; end

Matlab - randn

randn(n) : randomly produce an nxn matrix
using ... method of generator

randn(m,n) : randomly produce an mxn matrix
using ... method of generator

s = randn(method) returns in s the current internal state of the generator selected by method.
It does not change the generator being used.
'state' Uses Marsaglia's ziggurat algorithm (the default in MATLAB versions 5 and later).
The period is approximately 2^64.
'seed' Uses the polar algorithm (the default in MATLAB version 4).
The period is approximately (2^31-1)*(pi/8).
randn(method, n) initializes the state of this generator using the value of 100
Different initializers give different result.
Same initializer gives same result.

So, basically, if we want to keep the same "random" set of number through several calculation, we can do it by storing internal state.
Note that after each time we call randn, after the function generalize some random numbers, its internal state changes. If we want it to reproduce some set of "random" numbers again, we need to keep track of the internal state.

randn('state', 0) Set randn to its default initial state
randn('state', sum(100*clock)) Initialize randn to a different state each time

Example:
s = randn('state') % Save the current state
u1 = randn(3) % Generate 3x3 values
randn('state', s) % reset the state
u2 = randn(3) % contains exactly the same values as u1

Simulate SDEs using Matlab

http://www.maths.strath.ac.uk/~aas96106/algfiles.html

Numerical Computing with Matlab

Individual Chapters from Mathworks.