Last update 25.1.99

L2.html

3D-plotting in Matlab, meshgrid

Working with grids, interpolation

>> x=linspace(-2,2,6), y=linspace(-1,1,3) x = -2.0000 -1.2000 -0.4000 0.4000 1.2000 2.0000 y = -1 0 1
Consider f(x,y)=xy
>> [X,Y]=meshgrid(x,y) X = -2.0000 -1.2000 -0.4000 0.4000 1.2000 2.0000 -2.0000 -1.2000 -0.4000 0.4000 1.2000 2.0000 -2.0000 -1.2000 -0.4000 0.4000 1.2000 2.0000 Y = -1 -1 -1 -1 -1 -1 0 0 0 0 0 0 1 1 1 1 1 1 >> Z=X.*Y
Note: In a "real situation" remember the semicolon (;)
>> mesh(x,y,Z) >> surf(x,y,Z) >> hold on >> plot3(X(:),Y(:),Z(:),'o') % What do you think of this, nice trick?
In general, if you want to make a surface plot of z=f(x,y), you'd better write a function m-file. Like for instance if f(x,y)=exp(-(x-3)³ +(y-2)²), you edit a file:
%%%%%%%%%%%%%%% f.m %%%%%%%%%%%%%% function z=f(x,y) z=exp(-(x-3).^3+(y-2).^2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> y=linspace(0,4,20);x=linspace(0,.1,30); >> [X,Y]=meshgrid(x,y); >> Z=f(X,Y); >> mesh(x,y,Z) >> surf(x,y,Z) >> for i=linspace(0,360); view(i,60);i,pause; end % Rotate horizontally >> hold on >> plot3(X(:),Y(:),Z(:),'o') >> figure >> plot3(X(:),Y(:),Z(:),'o') >> grid >> plot3(X(:),Y(:),Z(:))

How dense grid is good ?
This leads us to
Interpolation

Polynomial and spline interpolation in 1-dim.

>> coeff=polyfit(xdata,ydata,n); x=linspace(a,b);y=polyval(coeff,x); >> ys=spline(xdata,ydata,xev);

Interpolating surface data

% Assume our function is known at some data points. (By for eg. FD or FEM- % calculation). % How do we evaluate it in intermediate points? [x,y]=meshgrid(-3:3); z=peaks(x,y); surf(x,y,z) [xi,yi]=meshgrid(-3:0.5:3); figure zi=interp2(x,y,z,xi,yi); % nearest, bilinear surf(xi,yi,zi) zi=interp2(x,y,z,xi,yi,'bicubic'); surf(xi,yi,zi) % We may also want to go in the reverse direction. We may have done a very % accurate computation, but there's no point using the mesh- or surf-functions % for a very dense data. (It takes very long and the result doesn't look good.) % Also, we may want to evaluate the error at certain data points that are a % subset of the computed values. % % Fine, it works both ways: zii=interp2(xi,yi,zi,x,y); surf(x,y,zii) % Our original "peaks".
See also:
DELAUNAY,VORONOI, TRIMESH, TRISURF, GRIDDATA, CONVHULL, DSEARCH, TSEARCH

Sparse matrices
Below is a list (a little shortened) of the command
help sparfun . Let's first look at:
SPDIAGS
It has 4 forms, let's have a look at 2 of them
1) A=spdiags(B,d,m,n) % B - matrix, columns are the diagonals % d - poiner vector [-1 0 1] would be tridiag. % m,n - size of matrix Example: N=5;yks=ones(N,1);ItoN=(1:N)';B1=[ItoN yks ItoN] B1 = 1 1 1 2 1 2 3 1 3 4 1 4 5 1 5 >> A=spdiags(B1,[-1 0 1],N,N);full(A) ans = 1 2 0 0 0 1 1 3 0 0 0 2 1 4 0 0 0 3 1 5 0 0 0 4 1 2) [B,d]=spdiags(A) B = 1 1 0 2 1 2 3 1 3 4 1 4 0 1 5 d = -1 0 1 Observations: ============= 1) The first n-1 elements of B(:,1) are used the last n-1 elements of B(:,3) are used. 2) Especially handy to create constant diagonals, no need to count lengths (just use the length of main diagonal for all subdiagonals as well) Sparse matrix functions. Elementary sparse matrices. speye - Sparse identity matrix. sprandn - Sparse random matrix. sprandsym - Sparse symmetric random matrix. spdiags - Sparse matrix formed from diagonals. <--- Full to sparse conversion. sparse - Create sparse matrix from nonzeros and indices. full - Convert sparse matrix to full matrix. find - Find indices of nonzero entries. spconvert - Convert from sparse matrix external format. Working with nonzero entries of sparse matrices. nnz - Number of nonzero entries. nonzeros - Nonzero entries. nzmax - Amount of storage allocated for nonzero entries. spones - Replace nonzero entries with ones. spalloc - Allocate memory for nonzero entries. issparse - True if matrix is sparse. spfun - Apply function to nonzero entries. Visualizing sparse matrices. spy - Visualize sparsity structure. gplot - Plot graph, as in "graph theory". Reordering algorithms. colmmd - Column minimum degree. symmmd - Symmetric minimum degree. symrcm - Reverse Cuthill-McKee ordering. colperm - Order columns based on nonzero count. randperm - Random permutation vector. dmperm - Dulmage-Mendelsohn decomposition. Norm, condition number, and rank. normest - Estimate 2-norm. condest - Estimate 1-norm condition. sprank - Structural rank. .....

Matlab's sparse matrix demos
Just give one of the following commands in Matlab
sparsity % The effect of reordering to fill-in % in Cholesky factorization buckydem % Connectivity graph sepdemo % Later? airfoil % Later?

How to solve a sparse system
One can use x=A\b but if the matrix has more structure (as is the case usually with FD or FEM-calculations), we should use something better.
Conjugate gradient for symmetric and positive definite

>> help pcg PCG Preconditioned Conjugate Gradients Method X = PCG(A,B) attempts to solve the system of linear equations A*X=B for X. The coefficient matrix A must be symmetric and positive definite and the right hand side (column) vector B must have length N, where A is N-by-N. PCG will start iterating from an initial guess which by default is an all zero vector of length N. Iterates are produced until the method either converges, fails, or has computed the maximum number of iterations. Convergence is achieved when an iterate X has relative residual NORM(B-A*X)/NORM(B) less than or equal to the tolerance of the method. The default tolerance is 1e-6. The default maximum number of iterations is the minimum of N and 20. No preconditioning is used. ...

ODE's in Matlab

help ode45
Give the Matlab-command helpdesk and type for instance ode45 in the search box.
See also here
To use the ODE-solvers you first need to be able to transform your equation(s) to a first order system (by hand). (If you have one 1^st order equation, then ther's no "paper/pencil work" of course.)
If your equation is y''+t*y'+2*y³=0 , you first write it as (Take y₁=y, y₂=y' .)
y₁'=y₂ y₂'=-2 y₁³-t y₂
In other words, the system is defined by
Y'=F(t,Y), where F(t,Y)=[Y(2);-2*Y(1)^3-t*Y(2)]
where we have used "semi Matlab syntax".
To define the system in Matlab we thus have to write the appropriate .m-file:
%%File F.m%%(You can't trust case sensitivity in file names,remember dos-files) z=F(t,Y) z=[Y(2);-2*Y(1)^3-t*Y(2)];
Note: This needn't be "array-smart" (it doesn't hurt if it is, though). The reason is that only scalar computations are used between the (scalar) components. (If you do implement Euler or Runge-Kutta or similar by yourself, you will immediately see this.)
Note: Even if your system were autonomous, you have to include t (or wahatever you want to call it) in the parameter list. Of course the parameter list must have a fixed form for the general ODEsolver to be able to call it.
Now call an ODE-solver
The solver can be ode45 or whichever ...
[T,Y]=solver('F',[0,5],[0;1]); time-span,IVals
The result is returned as a 3-col matrix: [T-col,Y1-col, Y2-col]. Plots can now be produced by ordinary plot as:
plot(T,Y) % time series plot(Y(:,1),Y(:,2)) % phase plane

Some Maple principles
Sources

[HAM] Heikki Apiola: Symbolista ja numeerista matematiikkaa Maple- ohjelmalla, Otatieto 588 Tässä on kirjaan liittyvä www-sivu.
Pikaopas
Quick reference Guide (Waterloo)

use maple xmaple&

Worksheet, using help (?)
expressions, assignment (:=), (semi)colon (;:), functions
sequences, lists, sets
Matrices and vectors: evalm(A&*B);
Solving diff. equ. (ODE and PDE) dsolve, dsolve(...numeric), DEtools
Graphics

Maple worksheet for separation of variables
Heat equ., download the worksheet

Open the worksheet in your Maple-session. It may be instructive to modify the worksheet into the case of Laplace equ, that's all you need for "Harj. 1" probl. 4.