>> 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 1Consider 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.*YNote: 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)3 +(y-2)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(:))
>> coeff=polyfit(xdata,ydata,n); x=linspace(a,b);y=polyval(coeff,x); >> ys=spline(xdata,ydata,xev);
% 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
help sparfun . Let's first look at:
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. .....
sparsity % The effect of reordering to fill-in % in Cholesky factorization buckydem % Connectivity graph sepdemo % Later? airfoil % Later?
>> 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. ...
If your equation is y''+t*y'+2*y3=0 , you first write it as (Take y1=y, y2=y' .)
y1'=y2 y2'=-2 y13-t y2In 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.
[T,Y]=solver('F',[0,5],[0;1]); time-span,IValsThe 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
use maple xmaple&
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.