Runge Kutta 现象

多项式插值一般在2-5次为宜，最多不超过6，7次。也就是说，在节点数很多的时候，最好将它们分组，分段，对每段则分别用低次多项式插值，这样效果比较接近真实。片面的追求高次多项式，将不会得到好结果。

一个著名的例子是Runge给的，他是Weierstruss在Berlin的学生，也是F. Klein在Gottingen的同事。
用于模拟常微分方程的解的重要的一类迭代法，龙格－库塔法（Runge-Kutta）是由他和Kutta共同发明的。月球上的Runge陨石坑 (Runge crater) 以他命名。

f(x)=1/(1+25*x^2)
在[-1，1]内取n+1（在这里取n=10）个等距节点，做n次牛顿插值的结果如下，

其中红线为插值结果，黑线为被插值的函数，可以看出，在[-.2，.2]区间内，插值结果很好，但随着区域的扩大，出现了很大的偏差，特别是在-1，1两点附近，插值多项式基本不能描述被插函数。
这就是Runge-Kutta现象。


源代码如下，
=========demo_RungeKutta.m:

x_node =-1:.2:1;
y_node = 1./(1+25.*x_node.*x_node);
a = mynewtonCoeff(x_node,y_node);
x=-1:.01:1;
y = mynewtonPoly(a,x_node,x);
plot(x_node,y_node,'*b',x,y,'r',x,y1,'k');

===========mynewtonCoeff.m:

function a = mynewtonCoeff(xData,yData)
% Returns coefficients of Newton's polynomial.
% USAGE: a = newtonCoeff(xData,yData)
% xData = x-coordinates of data points.
% yData = y-coordinates of data points.

n = length(xData);
a = yData;
for k = 2:n
a(k:n) = (a(k:n) - a(k-1))./(xData(k:n) - xData(k-1));
end

=========mynewtonPoly.m

function p = mynewtonPoly(a,xData,x)
% evaluate the value of Newton's polynomial at x.
% USAGE: p = newtonPoly(a,xData,x)
% a = coefficient array of the polynomial;
% must be computed first by newtonCoeff.
% xData = x-coordinates of data points.

n = length(xData);

p = a(n)*ones(size(x));
for k = 1:n-1;
p = a(n-k) + (x - xData(n-k)).*p;
end

Ax=b直接法

Why?
Ax=b是数值分析的基础，没有几个算法里面不需要解线性方程组的。

Linear, algebraic equations occur in almost all branches of numerical analysis. But their most visible application in engineering is in the analysis of linear systems (any
system whose response is proportional to the input is deemed to be linear). Linear systems include structures, elastic solids, heat flow, seepage of fluids, electromagnetic fields and electric circuits; i.e.,most topics taught in an engineering curriculum.

If the system is discrete, such as a truss（桁架（truss）是以剛性直桿為材料，在桿件之兩端以平滑之銷釘結合（pin joint）而成之剛性構架。在工程上運用甚廣，且具實用性及經濟性。桁架之桿件稱為構件（member），接合點亦稱為節點。實際上桁架之節點是用焊接而成，但為簡化設計，各直桿端部均假設以光滑銷釘連接，故各直桿均以二力桿件處理。） or an electric circuit, then its analysis leads directly to linear algebraic equations. In the case of a statically determinate truss, for example, the equations arise when the equilibrium conditions of the joints are written down. The unknowns x1, x2, . . . , xn represent the forces in the members and the support reactions, and the constants b1, b2, . . . , bn are the prescribed external
loads.

The behavior of continuous systems is described by differential equations, rather than algebraic equations. However, because numerical analysis can deal only with discrete variables, it is first necessary to approximate a differential equation with a
system of algebraic equations. The well-known finite difference, finite element and boundary element methods of analysis work in this manner. They use different approximations to achieve the “discretization,”but in each case the final task is the same:
solve a system(often a very large system) of linear, algebraic equations.

In summary, the modeling of linear systems invariably gives rise to equations of the form Ax = b, where b is the input and x represents the response of the system.

The coefficient matrix A, which reflects the characteristics of the system, is independent of the input. In other words, if the input is changed, the equations have to be solved again with a different b, but the same A. Therefore, it is desirable to have
an equation-solving algorithm that can handle any number of constant vectors with minimal computational effort.

解Ax=b的法子大概可分为两类，
1 直接法
2 迭代法

Overview of Direct Methods

Table 2.1 lists three popular direct methods, each of which uses elementary operations to produce its own final form of easy-to-solve equations.

Method Initial form Final form
Gauss elimination Ax = b Ux = c
LU decomposition Ax = b LUx = b
Gauss–Jordan elimination Ax = b Ix = c

LU法的优点是，对于某一个特定的A，当b变动时，仅需要forward substitution 和 back substitution就可以得到解。

function A = LUdec(A)
%Dolittle's decomposition
%returns A=[L\U]
%USAGE: A = LUdec(A)

n = size(A,1);
for i =1:n-1
for j =i+1:n
if(A(j,i))
l = A(j,i)/A(i,i);
A(j,i+1:n)= A(j,i+1:n)- l*A(i,i+1:n);
A(j,i)= l;
end
end
end

function x = LUsol(A,b)
% first call LUdec to decompose A=L*U
% then, Solves L*U*b = x by first forward substitution
% and finally back substitution
% USAGE: x = LUsol(A,b)

A = LUdec(A);

n=size(A,1);

for i = 2:n
b(i,:)= b(i,:) - A(i,1:i-1)*b(1:i-1,:);
end

for k = n:-1:1
b(k,:) = (b(k,:) - A(k,k+1:n)*b(k+1:n,:))/A(k,k);
end

x=b;

上面两个函数实现了LU分解和用LU的结果解方程。注意这里的代码可以解决Ax=b，当b不止一列的情形。

实现最简单的Gaussian Elimination

不考虑任何特殊情况，最简单的高斯消去法可实现如下，

% the most simple Gauss Elimination

function b = Gauss_Elim(A,b)

n = length(b);

%Elimination phase
for i =1:n-1
for j =i+1:n
if(A(j,i))
l = A(j,i)/A(i,i);
A(j,i:n)= A(j,i:n)- l*A(i,i:n);
b(j) = b(j)-l*b(i);
end
end
end

%back substitution
b(n)= b(n)/A(n,n);
for i = (n-1):-1:1
b(i)=b(i)-A(i,i+1:n)*b(i+1:n);
b(i)=b(i)/A(i,i);
end

试比较下面的代码，

function b = g(A,b)
n = length(b);
%-----------------Elimination phase-------------
for k = 1:n-1
for i = k+1:n
if A(i,k) ~= 0
lambda = A(i,k)/A(k,k);
A(i,k+1:n) = A(i,k+1:n) - lambda*A(k,k+1:n);
b(i)= b(i) - lambda*b(k);
end
end
end

size(A)

%--------------Back substitution phase-----------
for k = n:-1:1
k
b(k) = (b(k) - A(k,k+1:n)*b(k+1:n))/A(k,k)
end

这两段代码的区别，
1 在消去阶段，区别在于对主对角线元素下元素的处理。前者作了处理，后者未作处理。
2在回代阶段，后者代码中当k=n时，A(k,k+1:n)相当于A(n,n+1:n)，在matlab中实际上是一个Empty matrix: 1-by-0，
但是当A(k,k+1:n)*b(k+1:n)，结果会是0。这可能是matlab内部的一个设计，使得代码变得简洁。

sin(1/x)的内里乾坤

在学习微积分的时候，大多数教程都会给出 y=sin(1/x)，作为不连续函数的例子。这个函数的草图很容易画，大多数教程上也会有这么一个图，
x=-1:.0011:1;
y = sin(1./x);
plot(x,y,'-b');




从这图上看不出什么模式出来，然而，如果用散点图，
x=-1:.0011:1;
y = sin(1./x);
plot(x,y,'.r');
则如下图，

在散点图中似乎存在着某种规律和模式。

我们继续用下面的代码来仔细看看，

% this program show some interesting patterns hidden in the "randomness" of
% sin(1/x) near zero.

n=1:4000;
x = 1./n; % range of x is [.00025 1];
y = sin(n);
plot(x,y,'.b');
legend('x=1/n y=sin(n) n=1:4000','Location','NorthOutside');
legend boxoff;

这个图显然不太合理，因为大多数点都被压缩到很窄小的区域，而居中广阔的区域却只有少数几个点享用。因此，我们需要调整x轴的分辨率和定义范围。

如何调整？
我们可以用histgram函数来看看数据定义域的直方图，由此可以大概知道数据的分布情况。
用简单命令，
N =hist(x)

N =

3991 5 1 1 1 0 0 0 0 1

可以看出，x在[0,0.1]这个区间内，有3991个点！
而在[0.1,1.0]这个区间内只有9个点！
可见分布的极其不均匀。

因此，为了让“镜头”对准“大多数”，我们需要调整“镜头”的取向和焦距了。
这可以由命令
axis([xmin xmax ymin ymax]);
来调节。

我们先试着用
axis([0.0 .1 -1 1]);
得到的结果是这样的，

说明镜头对准的是少数的“权贵阶级”，仍未对准“沉默的大多数”，
实际上，用
xx= x(find(x<.1));hist(xx) 可以看出，在此区间内，分布极其不均匀， 此时的直方图如下，


















继续尝试，xx= x(find(x<.01));hist(xx) 发现此时分布有些好了，起码能看出点梯度来。


















此时的散点图如下，已能看出一些模式，但是大多数点仍偏居一隅！







再继续尝试，xx= x(find(x<.003));hist(xx)，

这时效果更为理想，有些接近正态分布了。 此时的散点图如下，





最后尝试，xx= x(find(x<.001));hist(xx)， 这时效果更为理想，更接近正态分布了。

此时的散点图如下，



所以我们最终决定采取x显示区间为[.0002， .001]，认为此时的效果最好。


从某种意义上，这是matlab的一个小小缺陷。因为axis auto并不能为我们选取一个适当的尺度来发现数据中隐藏的pattern，或许上面的分析能够改进这个功能的设计。

我们的观点是，合适的尺度应该让图中的主要点表现为类正态分布，如此，则人能够更好的把握数据中的pattern.

GeomView

GeomView and Ubuntu 7.10

---

Hi all,

I have installed GeomView on Ubuntu 7.10 with the command:

% sudo apt-get install geomview

Maxima

http://yenlung.math.nccu.edu.tw/linear/page1/page1.html

線性代數教材區

課程相關講義, 試題等連結或下載

Maixma在線性代數的應用 (PDF)

Maxima 的使用手冊, PDF 版

Maixma在線性代數的應用 (HTML)

Maxmia 的使用手冊, 線上版

1.1 Introduction

為什麼要學線性代數, 數學證明入門。

1.2 Vector Spaces

向量空間的基本定義及五個主要例子

1.3 Subspaces

子空間 (簡單的證明某個集合是向量空間)

1.4 Linear Combination and Systems of Linear Equations

線性組合和線性系統

1.5 Linear Dependent and Independent

線性相依和線性獨立

1.6 Bases and Dimension

基底是向量空間的主控者

2.1 Linear Transformations, Null Spaces, and Ranges

線性轉換和兩個重要的相關子空間

2.2 The Matrix Representation of a Linear Transformation

線性轉換用矩陣表示

2.3 Composition of Linear Transformation and Matrix Multiplication

兩個線性代數的合成, 對應的就是其表示矩陣相乘

2.4 Invertibility and Isomorphisms

基本列運算和基本矩陣, 計算的核心概念

3.1 Elementary Row Operations and Elementary Matrices

基本列運算和基本矩陣, 計算的核心概念

線性組合

Maxima 筆記系列

線性相依

Maxima 筆記系列

GNU TeXmacs Tutorial

http://www.texmacs.org/Tutorial/

Starting TeXmacs up

You may start TeXmacs by typing

    texmacs &

The first time that you start TeXmacs, the program will analyze your system in order to test whether all other packages which are needed by TeXmacs are present on your system. If no problems occur, then a window with the following aspect should appear on your screen:

Bit wise Operations

Bit wise Operations

Bitwise operations are operations which occur to the individual bits of a number. For example, 5 is written as 101 in base 2, that is, 5 = 101₂. Similarly, 3 = 11₂. In this case, operations occur on the individual bits.

For example, bit-wise and 3 and 5 is the and of the corresponding bits 101 and 11:

>> bitand(3, 5)

ans =

    1

>> bitor(3, 5)

ans =

    7

Relevant functions are:

• bitand
• bitor
• bitxor
• bitcmp
• bitshift
• bitset
• bitget
• bitmax

symbolic toolbox matlab

today install the symbolic toolbox.

=========

% here is a example to plot scalar field with surface and contour
% also the gradient and others
clear;
close all;

% another scalar field
[X,Y] = meshgrid(-15:.5:15, -10:.5:10);
Z = 5.0 * exp(-X.^2/90 - Y.^2/25);

%surface of Z(X,Y)
figure,surf(X,Y,Z),colorbar,axis equal

%contour
figure,
[c,h] = contour(X,Y,Z); clabel(c,h), colorbar,axis equal

%gradient
[px,py] = gradient(Z,.2,.15);
figure,
[c,h] = contour(X,Y,Z); colorbar,axis equal,hold on
%QUIVER(X,Y,U,V) plots velocity vectors as arrows with components (u,v)
% at the points (x,y).
quiver(X,Y,px,py), hold off, axis image

How to save image without border?

=======
Q:

Hi there,

Whenever I save a figure using:
imshow(I)
print frame1.tif -f1 -dtiff

or by using:
h = imshow(I)
saveas( h, 'frame1.tif', 'tif' )

, I always get a white border around the image (which is
present in the figure window). Does anyone know how to just
save the image itself with no border around it?

Thanks

===
A:

>>set(gca, 'position', [0 0 1 1], 'visible', 'off')
>> print(3,'-djpeg','C:\Documents and Settings\My Documents\3.jpg');


======
Also

用matlab保存图片的四种方法：

1 直接从菜单保存，有fig,eps,jpeg,gif,png,bmp等格式。

2 edit------〉copy figure，再粘贴到其他程序。

3 用saveas命令保存图片。

4 使用plot函数后紧接着用print函数。

print的三个参数：

（1）图形句柄，如果图形窗口标题栏是“Figure 3”,则句柄就是3.（2）单引号字符串，指定存储格式。

用gcf可以获取当前窗口句柄。

png格式：'-dpng'

jpeg: '-djpeg',

tiff: '-dtiff'

bmp: '-dbitmap'

（3）文件名。

详细用法请 help print

例：

>> x=-pi:2*pi/300:pi;
>> y=sin(x);
>> plot(x,y);
>> print(gcf,'-dpng','abc.png') % 保存为png格式的图片。

>> figure(2) % 新建一个句柄为2的图形窗口。
>> plot(x,cos(x)); % 在句柄为2的图形窗口上画图。
>> grid
>> print(2,'-djpeg','C:\abc.jpeg'); %将句柄为2的图形保存为jpeg/jpg格式的图片，文件名为'C:\abc.jpeg'。

FSPECIAL

help fspecial

FSPECIAL Create predefined 2-D filters.
H = FSPECIAL(TYPE) creates a two-dimensional filter H of the
specified type. Possible values for TYPE are:

'average' averaging filter
'disk' circular averaging filter
'gaussian' Gaussian lowpass filter
'laplacian' filter approximating the 2-D Laplacian operator
'log' Laplacian of Gaussian filter
'motion' motion filter
'prewitt' Prewitt horizontal edge-emphasizing filter
'sobel' Sobel horizontal edge-emphasizing filter
'unsharp' unsharp contrast enhancement filter

Depending on TYPE, FSPECIAL may take additional parameters
which you can supply. These parameters all have default
values.

H = FSPECIAL('average',HSIZE) returns an averaging filter H of size
HSIZE. HSIZE can be a vector specifying the number of rows and columns in
H or a scalar, in which case H is a square matrix.
The default HSIZE is [3 3].

H = FSPECIAL('disk',RADIUS) returns a circular averaging filter
(pillbox) within the square matrix of side 2*RADIUS+1.
The default RADIUS is 5.

H = FSPECIAL('gaussian',HSIZE,SIGMA) returns a rotationally
symmetric Gaussian lowpass filter of size HSIZE with standard
deviation SIGMA (positive). HSIZE can be a vector specifying the
number of rows and columns in H or a scalar, in which case H is a
square matrix.
The default HSIZE is [3 3], the default SIGMA is 0.5.

H = FSPECIAL('laplacian',ALPHA) returns a 3-by-3 filter
approximating the shape of the two-dimensional Laplacian
operator. The parameter ALPHA controls the shape of the
Laplacian and must be in the range 0.0 to 1.0.
The default ALPHA is 0.2.

H = FSPECIAL('log',HSIZE,SIGMA) returns a rotationally symmetric
Laplacian of Gaussian filter of size HSIZE with standard deviation
SIGMA (positive). HSIZE can be a vector specifying the number of rows
and columns in H or a scalar, in which case H is a square matrix.
The default HSIZE is [5 5], the default SIGMA is 0.5.

H = FSPECIAL('motion',LEN,THETA) returns a filter to approximate, once
convolved with an image, the linear motion of a camera by LEN pixels,
with an angle of THETA degrees in a counter-clockwise direction. The
filter becomes a vector for horizontal and vertical motions. The
default LEN is 9, the default THETA is 0, which corresponds to a
horizontal motion of 9 pixels.

H = FSPECIAL('prewitt') returns 3-by-3 filter that emphasizes
horizontal edges by approximating a vertical gradient. If you need to
emphasize vertical edges, transpose the filter H: H'.

[1 1 1;0 0 0;-1 -1 -1].

H = FSPECIAL('sobel') returns 3-by-3 filter that emphasizes
horizontal edges utilizing the smoothing effect by approximating a
vertical gradient. If you need to emphasize vertical edges, transpose
the filter H: H'.

[1 2 1;0 0 0;-1 -2 -1].

H = FSPECIAL('unsharp',ALPHA) returns a 3-by-3 unsharp contrast
enhancement filter. FSPECIAL creates the unsharp filter from the
negative of the Laplacian filter with parameter ALPHA. ALPHA controls
the shape of the Laplacian and must be in the range 0.0 to 1.0.
The default ALPHA is 0.2.

Class Support
-------------
H is of class double.

Example
-------
I = imread('cameraman.tif');
subplot(2,2,1);imshow(I);title('Original Image');
H = fspecial('motion',20,45);
MotionBlur = imfilter(I,H,'replicate');
subplot(2,2,2);imshow(MotionBlur);title('Motion Blurred Image');
H = fspecial('disk',10);
blurred = imfilter(I,H,'replicate');
subplot(2,2,3);imshow(blurred);title('Blurred Image');
H = fspecial('unsharp');
sharpened = imfilter(I,H,'replicate');
subplot(2,2,4);imshow(sharpened);title('Sharpened Image');

See also conv2, edge, filter2, fsamp2, fwind1, fwind2, imfilter.

Reference page in Help browser
doc fspecial

gradient

code

[x,y] = meshgrid(-2:.2:2,-1:.15:1);
z = x .* exp(-x.^2 - y.^2); [px,py] = gradient(z,.2,.15);
contour(x,y,z), hold on
quiver(x,y,px,py), hold off, axis image





------

Geometric Spatial Transformation

[x y 1]=[w z 1]T

T = [ t11 t12 0]
[ t21 t22 0]
[ t31 t32 1]

relational operators

关系运算

>> A

A =

1 2
3 4

>> B

B =

1 3
4 2

>> A < B

ans =

0 1
1 0



>> A >= 2

ans =

0 1
1 1

index 问题

对于二维矩阵，可以用二维的下标引用，这是最直观，自然的。但是在内存中，matlab矩阵的存储次序如何呢，是一列一列的存的。比如，

>> A = rand(2,5)

A =

0.4089 0.1436 0.0832 0.0304 0.7000
0.7080 0.8713 0.4617 0.7532 0.2145

在内存中，则是，

起始地址 Addr

Addr ：0.4089
Addr +1：0.7080
Addr +2：0.1436
Addr +3：0.8713
Addr +4：0.0832
Addr +5：0.4617
。。。

所以，
>> A(6)

ans =

0.4617

>> A(2,3)

ans =

0.4617

MATLAB functions

MATLAB provides a large number of standard elementary mathematical functions, including abs, sqrt, exp, and sin. Taking the square root or logarithm of a negative number is not an error; the appropriate complex result is produced automatically. MATLAB also provides many more advanced mathematical functions, including Bessel and gamma functions. Most of these functions accept complex arguments. For a list of the elementary mathematical functions, type

help elfun

For a list of more advanced mathematical and matrix functions, type

help specfun
help elmat

Some of the functions, like sqrt and sin, are built in. Built-in functions are part of the MATLAB core so they are very efficient, but the computational details are not readily accessible. Other functions, like gamma and sinh, are implemented in M-files.

There are some differences between built-in functions and other functions. For example, for built-in functions, you cannot see the code. For other functions, you can see the code and even modify it if you want.

Several special functions provide values of useful constants.

pi


3.14159265...

i


Imaginary unit,

j


Same as i

eps


Floating-point relative precision,

realmin


Smallest floating-point number,

realmax


Largest floating-point number,

Inf


Infinity

NaN


Not-a-number

Infinity is generated by dividing a nonzero value by zero, or by evaluating well defined mathematical expressions that overflow, i.e., exceed realmax. Not-a-number is generated by trying to evaluate expressions like 0/0 or Inf-Inf that do not have well defined mathematical values.

The function names are not reserved. It is possible to overwrite any of them with a new variable, such as

eps = 1.e-6

and then use that value in subsequent calculations. The original function can be restored with

clear eps


>> help elfun
Elementary math functions.

Trigonometric.
sin - Sine.
sind - Sine of argument in degrees.
sinh - Hyperbolic sine.
asin - Inverse sine.
asind - Inverse sine, result in degrees.
asinh - Inverse hyperbolic sine.
cos - Cosine.
cosd - Cosine of argument in degrees.
cosh - Hyperbolic cosine.
acos - Inverse cosine.
acosd - Inverse cosine, result in degrees.
acosh - Inverse hyperbolic cosine.
tan - Tangent.
tand - Tangent of argument in degrees.
tanh - Hyperbolic tangent.
atan - Inverse tangent.
atand - Inverse tangent, result in degrees.
atan2 - Four quadrant inverse tangent.
atanh - Inverse hyperbolic tangent.
sec - Secant.
secd - Secant of argument in degrees.
sech - Hyperbolic secant.
asec - Inverse secant.
asecd - Inverse secant, result in degrees.
asech - Inverse hyperbolic secant.
csc - Cosecant.
cscd - Cosecant of argument in degrees.
csch - Hyperbolic cosecant.
acsc - Inverse cosecant.
acscd - Inverse cosecant, result in degrees.
acsch - Inverse hyperbolic cosecant.
cot - Cotangent.
cotd - Cotangent of argument in degrees.
coth - Hyperbolic cotangent.
acot - Inverse cotangent.
acotd - Inverse cotangent, result in degrees.
acoth - Inverse hyperbolic cotangent.
hypot - Square root of sum of squares.

Exponential.
exp - Exponential.
expm1 - Compute exp(x)-1 accurately.
log - Natural logarithm.
log1p - Compute log(1+x) accurately.
log10 - Common (base 10) logarithm.
log2 - Base 2 logarithm and dissect floating point number.
pow2 - Base 2 power and scale floating point number.
realpow - Power that will error out on complex result.
reallog - Natural logarithm of real number.
realsqrt - Square root of number greater than or equal to zero.
sqrt - Square root.
nthroot - Real n-th root of real numbers.
nextpow2 - Next higher power of 2.

Complex.
abs - Absolute value.
angle - Phase angle.
complex - Construct complex data from real and imaginary parts.
conj - Complex conjugate.
imag - Complex imaginary part.
real - Complex real part.
unwrap - Unwrap phase angle.
isreal - True for real array.
cplxpair - Sort numbers into complex conjugate pairs.

Rounding and remainder.
fix - Round towards zero.
floor - Round towards minus infinity.
ceil - Round towards plus infinity.
round - Round towards nearest integer.
mod - Modulus (signed remainder after division).
rem - Remainder after division.
sign - Signum.

>> help elmat
Elementary matrices and matrix manipulation.

Elementary matrices.
zeros - Zeros array.
ones - Ones array.
eye - Identity matrix.
repmat - Replicate and tile array.
rand - Uniformly distributed random numbers.
randn - Normally distributed random numbers.
linspace - Linearly spaced vector.
logspace - Logarithmically spaced vector.
freqspace - Frequency spacing for frequency response.
meshgrid - X and Y arrays for 3-D plots.
accumarray - Construct an array with accumulation.
: - Regularly spaced vector and index into matrix.

Basic array information.
size - Size of array.
length - Length of vector.
ndims - Number of dimensions.
numel - Number of elements.
disp - Display matrix or text.
isempty - True for empty array.
isequal - True if arrays are numerically equal.
isequalwithequalnans - True if arrays are numerically equal.

Matrix manipulation.
cat - Concatenate arrays.
reshape - Change size.
diag - Diagonal matrices and diagonals of matrix.
blkdiag - Block diagonal concatenation.
tril - Extract lower triangular part.
triu - Extract upper triangular part.
fliplr - Flip matrix in left/right direction.
flipud - Flip matrix in up/down direction.
flipdim - Flip matrix along specified dimension.
rot90 - Rotate matrix 90 degrees.
: - Regularly spaced vector and index into matrix.
find - Find indices of nonzero elements.
end - Last index.
sub2ind - Linear index from multiple subscripts.
ind2sub - Multiple subscripts from linear index.
bsxfun - Binary singleton expansion function.

Multi-dimensional array functions.
ndgrid - Generate arrays for N-D functions and interpolation.
permute - Permute array dimensions.
ipermute - Inverse permute array dimensions.
shiftdim - Shift dimensions.
circshift - Shift array circularly.
squeeze - Remove singleton dimensions.

Array utility functions.
isscalar - True for scalar.
isvector - True for vector.

Special variables and constants.
ans - Most recent answer.
eps - Floating point relative accuracy.
realmax - Largest positive floating point number.
realmin - Smallest positive floating point number.
pi - 3.1415926535897....
i - Imaginary unit.
inf - Infinity.
nan - Not-a-Number.
isnan - True for Not-a-Number.
isinf - True for infinite elements.
isfinite - True for finite elements.
j - Imaginary unit.
why - Succinct answer.

Specialized matrices.
compan - Companion matrix.
gallery - Higham test matrices.
hadamard - Hadamard matrix.
hankel - Hankel matrix.
hilb - Hilbert matrix.
invhilb - Inverse Hilbert matrix.
magic - Magic square.
pascal - Pascal matrix.
rosser - Classic symmetric eigenvalue test problem.
toeplitz - Toeplitz matrix.
vander - Vandermonde matrix.
wilkinson - Wilkinson's eigenvalue test matrix.

>> help specfun
Specialized math functions.

Specialized math functions.
airy - Airy functions.
besselj - Bessel function of the first kind.
bessely - Bessel function of the second kind.
besselh - Bessel functions of the third kind (Hankel function).
besseli - Modified Bessel function of the first kind.
besselk - Modified Bessel function of the second kind.
beta - Beta function.
betainc - Incomplete beta function.
betaln - Logarithm of beta function.
ellipj - Jacobi elliptic functions.
ellipke - Complete elliptic integral.
erf - Error function.
erfc - Complementary error function.
erfcx - Scaled complementary error function.
erfinv - Inverse error function.
expint - Exponential integral function.
gamma - Gamma function.
gammainc - Incomplete gamma function.
gammaln - Logarithm of gamma function.
psi - Psi (polygamma) function.
legendre - Associated Legendre function.
cross - Vector cross product.
dot - Vector dot product.

Number theoretic functions.
factor - Prime factors.
isprime - True for prime numbers.
primes - Generate list of prime numbers.
gcd - Greatest common divisor.
lcm - Least common multiple.
rat - Rational approximation.
rats - Rational output.
perms - All possible permutations.
nchoosek - All combinations of N elements taken K at a time.
factorial - Factorial function.

Coordinate transforms.
cart2sph - Transform Cartesian to spherical coordinates.
cart2pol - Transform Cartesian to polar coordinates.
pol2cart - Transform polar to Cartesian coordinates.
sph2cart - Transform spherical to Cartesian coordinates.
hsv2rgb - Convert hue-saturation-value colors to red-green-blue.
rgb2hsv - Convert red-green-blue colors to hue-saturation-value.

Save workspace variables to disk

save -
Save workspace variables to disk

Graphical Interface

As an alternative to the save function, select Save Workspace As from the File menu in the MATLAB® desktop, or use the Workspace browser.

Syntax

save
save filename
save filename content
save filename options
save filename content options
save('filename', 'var1', 'var2', ...)


Examples

Example 1
Save all variables from the workspace in binary MAT-file test.mat:
save test.mat
Example 2
Save variables p and q in binary MAT-file test.mat.
In this example, the file name is stored in a variable, savefile. You must call save using the function syntax of the command if you intend to reference the file name through a variable.
savefile = 'test.mat';
p = rand(1, 10);
q = ones(10);
save(savefile, 'p', 'q')
Example 3
Save the variables vol and temp in ASCII format to a file named june10:
save('d:\mymfiles\june10','vol','temp','-ASCII')
Example 4
Save the fields of structure s1 as individual variables rather than as an entire structure.
s1.a = 12.7; s1.b = {'abc', [4 5; 6 7]}; s1.c = 'Hello!';
save newstruct.mat -struct s1;
clear
Check what was saved to newstruct.mat:
whos -file newstruct.mat
Name Size Bytes Class

a 1x1 8 double array
b 1x2 158 cell array
c 1x6 12 char array

Grand total is 16 elements using 178 bytes
Read only the b field into the MATLAB workspace.
str = load('newstruct.mat', 'b')
str =
b: {'abc' [2x2 double]}
Example 5
Using regular expressions, save in MAT-file mydata.mat those variables with names that begin with Mon, Tue, or Wed:
save('mydata', '-regexp', '^Mon|^Tue|^Wed');
Here is another way of doing the same thing. In this case, there are three separate expression arguments:
save('mydata', '-regexp', '^Mon', '^Tue', '^Wed');
Example 6
Save a 3000-by-3000 matrix uncompressed to file c1.mat, and compressed to file c2.mat. The compressed file uses about one quarter the disk space required to store the uncompressed data:
x = ones(3000);
y = uint32(rand(3000) * 100);

save -v6 c1 x y % Save without compression
save -v7 c2 x y % Save with compression

d1 = dir('c1.mat');
d2 = dir('c2.mat');

d1.bytes
ans =
45000240 % Size of the uncompressed data in bytes.
d2.bytes
ans =
11985283 % Size of the compressed data in bytes.

d2.bytes/d1.bytes
ans =
0.2663 % Ratio of compressed to uncompressed

tips

You can also assign a scalar to all entries of a submatrix. Try:
A(:, [2 4]) = 99

You can delete rows or columns of a matrix by assigning the empty matrix ([]) to them. Try:
A(:, [2 4]) = []


>> x= rand(1,5)

x =

0.7513 0.2551 0.5060 0.6991 0.8909

>> x = x(end:-1:1)

x =

0.8909 0.6991 0.5060 0.2551 0.7513

Basic 2-d Plot Commands

Basic 2-d Plot Commands

The basic command for plotting in Matlab is: plot(x,y). The important thing to remember is that in Matlab, variables such as x and y are matrices, which must be of the same size. The plot command simply creates a figure window and plots the data in the variable x versus those in y and interpolates linearly between the points. If x is a vector and y are a matrix, then plot(x,y) will plot the rows or columns of y versus the vector x, depending on whether the row or column length matches the length of x.

We are interested in using Matlab to plot functions. The important point is: all built in functions in Matlab are vectorized. This means that if x is a matrix, then the command sin(x) will create a matrix of the same size whose entries are computed by taking the sine of each entry of x. Hence, to create a plot of sin(x) on say the interval [-pi,pi] we would use:

>>x=-pi:pi/20:pi;
>>plot(x,sin(x))

The first line creates a vector of data for the independent variable; the step size is pi/20. To find out how long the vector is type

>>length(x) % no semicolon produces output

ans = 41 % this is the form of output

Alternatively, we could type

>>size(x)

ans = 1 41 % x is a row vector of length 41

Multiple plots on the same axis can be done in two ways.

>>x=-pi:pi/20:pi;
>>plot(x,sin(x))
>>hold on % this command forces the same set of axis for further plots
>>plot(x,cos(x))
>>grid on % places a grid on the axis
>>hold off % future plot commands will go to a new set of axis

Alternatively, the same thing can be accomplished with the command

>>plot(x,sin(x),x,cos(x))

Either way produces something like the following output.



Alternatively, one may wish to plot the two functions in the same figure window but on a different set of axes. The command is

>>subplot(1,2,1),plot(x,sin(x))
>>grid on
>>subplot(1,2,2),plot(x,cos(x))
>>grid on

NOTE: If you forget to type grid on, or any other axis command below after the first plot, simply click on the axis region and then type the command.

The output looks like:



Several observations are in order:

The x-axis ranges from -4 to 4, not matching the original data.
Except for color, we cannot distinguish between the graphs.
The axes are not labeled.
Are there other line styles?
There is no title on the figure.
Can we control the tick mark positions on the x and y axes?
Axis commands:

axis tight (Sets axis limits to the range of the plotted data).
axis([xmin xmax ymin ymax]) (Used before the plot command, sets the limits as given in the vector)
axis equal (Set the aspect ratio so that equal tick mark increments have the same size on each axis)
axis normal (Restore the axis box to full size and removes any restriction on scaling)
axis off (Turn off the axis, tick marks, and background)
axis on (Turn on the axis, tick marks, and background)
Line styles. The curves depicted in a given set of axis can have a variety of line styles. For example

>>plot(x,sin(x),x,cos(x),'--') % plots sin(x) and cos(x) on the same axis with cosine plotted with a dashed line

>>plot(x,sin(x),'LineWidth',2) % plots sin(x) with line thickness of 2 points; default is 0.5; possible choices 2,3,4.

Labeling the axis:

xlabel('string') % string represents the desired label
ylabel('string')
zlabel('string') % For use in 3-d plots
Annotating the plot:

legend('string1','string2',.....,m)
The strings string1, string2, etc are the desired labels for the different functions plotted. The variable m placed at the end specifies the position of the legend box: m=1, upper right hand corner, m=2 upper left hand corner, m=3 lower left hand corner, m=4 lower right hand corner, m=-1 to the right of the axis, without m specified the default is m=1.
gtext('string')
This command transforms the mouse pointer into crosshairs over the plot, when clicked over a desired position, the string will appear on the plot.
Plot title command:

title('string') % places string above the plot
gtext can be used as above for manual positioning of a title.
Tick mark values can be set with the command:

set(gca,'xtick',[vector of values])
gca means : get current axis object
'xtick' : could be 'ytick' or 'ztick'
[vector of values] : you specify the numerical values for the ticks.

Example: Sine & Cosine Plots

This example makes use of a variety of the above commands.

>>x=-pi:pi/20:pi;
>>plot(x,sin(x),x,cos(x),'--',2)
>>axis tight
>>legend('sin(x)','cos(x)',2)
>>title('An Example Plot')
>>xlabel('x'); ylabel('y')
>>set(gca,'xtick',[-3 0 3])
>>grid on

The output should look something like:



NOTE: Many of the plot refinements can be done from the figure window. Look under the Tool menu in the figure window. The tool bar also offers a variety of ways to edit the plot.

Example: Roots of an equation

Here we will obtain a plot of tan(pi*x) and -x. This is a graphical depiction of the roots of the equation tan(pi*x)=-x.

» x=-0.495:0.005:0.495; % vector for x
» x=[x x+1 x+2 x+3 x+4]; % making use of periodicity and avoiding singularities
» plot(x,tan(pi*x),'LineWidth',2)
» hold on
» axis([0 4.5 -6 4]) % adjusting the axis limits
» plot(x,-x) % plotting the line
» plot(x,zeros(length(x))) % plotting a horizontal line
» title('tan\pi\lambda = -\lambda')
» title('tan\pi\lambda = -\lambda : Graphical Depiction of Roots')
» gtext('\lambda_{1}')
» gtext('\lambda_{2}')
» gtext('\lambda_{3}')
» gtext('\lambda_{4}')

NOTE: In the last five lines we have used a LaTex interpreter built into Matlab for writing Greek letters. The backslash \ invokes the interpreter, LaTex command for Greek letters is to spell them, and the _{k} forces a subscript. The output looks as follows.



WARNING: In Matlab, the operand +, -, *, / are all matrix operations.

There is no problem with the first two: if x is a matrix, then x+5 simply adds 5 to every entry in x. However, if x and y are matrices of the same size, x*y performs matrix multiplication. If you want multiplication to be done element by element, the command is x.*y.

For example, if you want to plot the polynomial function (x^4-1)(x+3) on the interval [-4,2] one could write

>>x=-4:0.1:2;
>>plot(x,(x.^4-1).*(x+3))

验证一个命题

假设，C为mxn矩阵
则有
trace(C'*C)=norm(C,'fro');


------


>> C = rand(3,3)

C =

0.4898 0.7094 0.6797
0.4456 0.7547 0.6551
0.6463 0.2760 0.1626

>> trace(C'*C)

ans =

2.9227

>> norm(C,'fro')*norm(C,'fro')

ans =

2.9227


--------


help norm
NORM Matrix or vector norm.
For matrices...
NORM(X) is the largest singular value of X, max(svd(X)).
NORM(X,2) is the same as NORM(X).
NORM(X,1) is the 1-norm of X, the largest column sum,
= max(sum(abs(X))).
NORM(X,inf) is the infinity norm of X, the largest row sum,
= max(sum(abs(X'))).
NORM(X,'fro') is the Frobenius norm, sqrt(sum(diag(X'*X))).
NORM(X,P) is available for matrix X only if P is 1, 2, inf or 'fro'.

For vectors...
NORM(V,P) = sum(abs(V).^P)^(1/P).
NORM(V) = norm(V,2).
NORM(V,inf) = max(abs(V)).
NORM(V,-inf) = min(abs(V)).

See also cond, rcond, condest, normest, hypot.

Overloaded methods:
lti/norm
distributed/norm
mfilt.norm
adaptfilt.norm
idmodel/norm
dfilt.norm

Reference page in Help browser
doc norm