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不止一列的情形。
Sunday, October 11, 2009
实现最简单的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内部的一个设计,使得代码变得简洁。
% 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内部的一个设计,使得代码变得简洁。
Saturday, October 3, 2009
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.


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.


Monday, July 6, 2009
GeomView
GeomView and Ubuntu 7.10
Hi all,
I have installed GeomView on Ubuntu 7.10 with the command:
% sudo apt-get install geomview
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
線性代數教材區
課程相關講義, 試題等連結或下載
Maxima 的使用手冊, PDF 版
Maxmia 的使用手冊, 線上版
為什麼要學線性代數, 數學證明入門。
向量空間的基本定義及五個主要例子
子空間 (簡單的證明某個集合是向量空間)
線性相依和線性獨立
基底是向量空間的主控者
線性轉換和兩個重要的相關子空間
兩個線性代數的合成, 對應的就是其表示矩陣相乘
基本列運算和基本矩陣, 計算的核心概念
基本列運算和基本矩陣, 計算的核心概念
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:
Monday, June 15, 2009
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 = 1012. Similarly, 3 = 112. 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
Subscribe to:
Posts (Atom)