MATLAB-数值方法

Matlab_5 数值方法
多项式
线性方程组
代数方程
数据统计
数据插值
数据拟合
数值积分
数值微分
常微分方程
偏微分方程
多项式
Matlab提供了一组函数用于处理多项式运算。
（不适合处理大于10的高阶多项式）
常用的函数：roots, polyval, polyfit, …
%%多项式表达：
% 用一个行向量表示各阶系数，阶次降序排列。
p1 = [1 -3 2] % x2 - 3x + 2
p2 = [1 0 -2 3] % x3 - 2x + 3
%% 多项式求值：
polyval( p1, 2 )
polyval( p2, 1 )
% 绘制多项式x3 - 2x + 3的图形
x = -2:0.1:2;
y = polyval( p2, x);
plot(x,y,’.-’)
%% 多项式求根：
roots(p1)
roots(p2)
roots( [1 1 0 0] )
%% 构建多项式：
p = 1:5
r = roots( p )
pp = poly( r )
pp - p
% 由于截断误差，函数poly生成的系数有微小的偏差
% 有时候结果出现复数，可以用real函数提取实部，消除虚部的影响。
% 多项式加减运算：
%多项式加减就是系数向量的加减。
% 多项式乘法：（系数向量的卷积运算）
a = [ 1 3 5 7 ]
p = conv( a, 1:3 )
% 多项式除法：
[ q, r ] = deconv( p, a )
% q商, r余数
[ q, r ] = deconv( 2:6, a)
conv(a, q) + r
% 多项式微分：
a = 1:3, b = [1 1 1 1 1]
polyder( a )
polyder( b )
polyder( a, b) % 多项式乘积的导数。
polyder( conv(a,b) )
% 有理多项式<b/a>的导数，返回分子和分母多项式。
% 注意Matlab函数的输入输出参数
[ bb, aa ] = polyder( b, a)
[ bb, aa ] = polyder( [1 0], [1 -1] )
% 多项式积分：
p = polyint( [ 1 1 1 ] )
polyder( p )
polyint( [ 1 1 1 ], 10 ) % 指定积分常数为10
%%多项式拟合
%格式：[p, s, mu] = polyfit( x, y, n )
%% 例：
x = 0:0.1:1 ;
y = [ -0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2 ];
plot( x, y, ‘or’ )
n = 2;
p = polyfit(x, y, n) % 二阶拟合
x2 = linspace(0, 1);
y2 = polyval(p, x2);
hold on, plot( x2, y2, ‘.-’ ), hold off
p10 = polyfit(x, y, 10); % 十阶拟合
p10.’
y10 = polyval(p10, x2);
hold on, plot(x2, y10, ‘-m’ ), hold off
%%查看拟合效果
[p, s] = polyfit(x, y, 2)
yp = polyval(p, x);
sqrt( sum( (yp-y).^2 ) ) % s.normr
R2 = 1 - sum( (yp-y).^2 )/sum( y.^2 )
%%尺度变换的调用格式
[p, s, mu] = polyfit( x, y, 2)
mean(x) % mu(1)
std(x) % mu(2)
xp = (x-mu(1))/mu(2);
yp = polyval(p, xp);
plot( xp, y, ‘o’, xp, yp, ‘.-’ );
线性方程组
一般形式：Ax=B（B可以是矩阵）
齐次性：系数矩阵B是否为零向量
相容性：系数矩阵A与增广矩阵[A,B]的秩相等
满秩性：A的秩与未知数相等
% [m, n] = size(A); % m方程数目，n未知数数目
% freedom = n - rank(A) % 基本解系数目

norm(B) == 0
N非齐次方程 Y齐次方程
rank(A) == rank([A,B]) rank(A) == n
Y相容（有解） N不相容 Y零解 N无穷多解
rank(A) == n
最小二乘解
x = A\B
=pinv(A)B

（超定）

x = zeros(n,1)

x = null(A)
Y唯一解 N无穷多解: xt + x*
x = A\B
或
x = inv(A).B
（适定） 特解：xt = A\B
或xt = pinv(A).B
通解：x* = null(A)
（欠定）
欠定问题：pinv给出最小范数特解，A\B最多零元素特解。
方程组数目m不能决定方程解的特性，因为可能包含多余的方程式。
%例1
a = magic(8);
b = 260*ones(8,1);
a (:, 6:8) = [ ]
[n_equ, n_x] = size(a)
norm(b)
rank(a)
rank([a, b])
x0 = null(a)
xt1 = pinv(a)*b
xt2 = a\b
%例2
a = [1 0; 1 1; 2 1], b = (1:3).’
rank(a), rank([a b])

null(a)
inv(a)b
a\b
pinv(a)b
%% 矩阵特征值和特征向量
% 特征多项式：det( A - lambdaE)=0，特征值lambda
a = magic(3)
lambda = eig(a) %特征值
[v d] = eig(a) %特征向量和特征值对角阵av = v*d
poly(a) %特征多项式系数
roots(poly(a)) %特征值
%% 向量V的范数 ||V||n
% norm(V, n) %= ( sum(abs(V).^n) )^(1/n) % n缺省为2-范数
%%矩阵A的范数
% norm(A) % 2-范数： ||A||2 = sqrt(max(lambda(A’.A)))
%列范数：每列绝对值之和的最大值；行范数…
%% 逆阵和伪逆阵
% 逆阵：AB=BA=E; det(A) ~= 0
% 伪逆阵：ABA = A, BAB = B; rank(A) < n
% 矩阵的条件数：cond(A) = norm(A)*norm(A^-1)
%% Cram法则: Ax = B
% Ai = A; Ai(:, i) = B;
% xi = Di/D, D = |A|, Di = |Ai|
% 比较Cram法则的求解效率
n = 10;
a = rand(n); b = rand(n, 1);
a1 = a;
a1(:,1) = b;
tic, x_1 = det(a1)./det(a); toc %求解x1
tic, x_all = a\b; toc %求解全部x
超定问题：最小二乘法

%%将超定问题转换为适定问题求解
% 超定->适定（m行数，方程数目；n列数，未知数目）

a11x1 +… a1nxn = b1
…
am1x1 + … amnxn = bn

二范数最小化
ϵ^2=∑_(i=1)m▒∑_(j=1)^n▒(a_ij x_j-b_i )^2
〖∂/(∂x_k ) ϵ〗^2=0; k=1,2,…,n

〖∂/(∂x_k ) ϵ〗^2=∑_(i=1)m▒[2a_ik ∑_(j=1)^n▒〖(a_ij x_j-b_i)〗]
(k = 1, 2 … n)

%%例：a.x = b -> aa.x = bb

rng(0); a = rand(6,2); b = rand(6,1);

[m, n]=size(a); aa=zeros(n); bb=zeros(n,1);

for k=1:n
for j=1:n
aa(k,j)=0; %?
for i=1:m
aa(k,j)=aa(k,j)+a(i,k)*a(i,j);
end;
end;
bb(k)=0; %?
for i=1:m
bb(k)=bb(k)+a(i,k)*b(i);
end;
end;

x1 = aa\bb
x0 = a\b
代数方程（优化）
函数y = f(x)在x取何值时，能够得到特定的y值或极值，Matlab称为优化。
对于实际问题直接利用反函数f -1(y)确定x是很困难的，通常需要迭代求解。
这个迭代过程实际上是求解y – f(x) = 0的过程。
% humps函数
humps(x)=1/((x-0.2)^{2+0.01)+1/((x-0.9)}2+0.04)-6

x = linspace( -0.5, 1.5);
y = humps(x);
plot( x, y, ‘.-’);
grid on;
一元函数寻零—fzero
% 指定点附近搜索
fzero(@humps, 1) % 在x = 1附近开始搜索
[x, r, exitflag, output] = fzero(@humps, -1) % 返回值r为偏差
% 指定搜索区域
fzero( @humps, [1 2] ) % 区域两段函数值符号相反
fzero( @humps, [0 1] )
fzero( @humps, [-3 3] )
% 待处理的函数需在搜索区域连续
fzero( @tan, 0.1)
fzero( @tan, [0.1, 3] )
[x, r] = fzero( @tan, [0.1, 3] )
求解选项的设置和获取—optimset和optimget函数
help optimset % 查看帮助文件
%
options = optimset(‘Display’, ‘iter’, ‘TolX’, 1e-3)
[x, value] = fzero(@humps, [-1 2]