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], options )
一维最小值—fminbnd
调用格式：
[ xmin, value, exitflag, output ] = fminbnd( fun, x1, x2, options )
% 指定在[x1, x2]区间内搜索最小值。
% 在x = 0.5~0.8区域搜索humps函数最小值
options = optimset(‘Display’, ‘iter’);
[ xmin, value, flag, out ] = fminbnd( @humps, 0.5, 0.8, options )
% 恰当选择搜索区域
[ xmin, value, flag ] = fminbnd( @humps, 0.5, 1.4, options )
% 最大值问题
% 等同于寻找函数–f(x)的最小值。
h = @(x) -1humps(x) ;
[ xmax, value, flag ] = fminbnd( h, 0.5, 1.4, options )
多维最小值—fminsearch函数
% Rosenbrock函数和图形
f(x,y)=100(y-x^2 )^{2+〖(1-x)〗}2
[x, y] = meshgrid( -1.5:0.1:1.5, 0:0.1:2.0);
z = 100(y-x.^2).2 + (1-x).^2;
surf(x,y,z);
% 函数fminsearch输入函数的自变量写成向量形式
% 对于较复杂的待求解函数，最好写成M函数
h = @(x)( 100*(x(2)-x(1).^2).2 + (1-x(1)).^2 ) ;
[ xymin, value, flag, output ] = fminsearch( h, [-1, 2] ) % 在(-1,2)附近搜索
% 非线性拟合：fminsearch+norm
% 例：求解线性方程（仅作说明用）
% 用fminsearch最小化norm(Ax-B)，可求出xi，即最小二乘解
rng(0);
a = rand(5,3);
b = rand(5,1);
h = @(x) norm( a*x - b );
x = fminsearch(h, [1 1 1].’)
% 对数据系列(xi,yi)，由公式y=f(x,a)进行非线性拟合
%用函数fminsearch处理2范数norm(yfit-ydata)，得到系数a
%% 例：用公式拟合 y = a1exp(a2x)下列数据
x = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3];
y = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
plot(x, y, ‘or’)

h = @(a, x) a(1).exp(a(2).x);
h1 = @(a)norm(h(a,x)-y);
a0 = [100 0];
a = fminsearch(h1, a0)
xx = 0:90;
yy = h(a, xx);
hold on; plot(xx, yy); hold off
% …
求解非线性代数方程组—fsolve
doc fsolve % 查看函数fsolve用法
% 调用格式：[ x, fval, exitflag, output, jacobian] = fsolve( fun, x0, options )
% 待求解问题必须写成标准形式：F(x) = 0
% 例：求解方程组
2x1 – x2 = e-x1
-x1 + 2x2 = e-x2
% 例 ---------------------------------------------------------
function F = funfsolve( x )
F = [ 2x(1) - x(2) - exp(-x(1))
-x(1) + 2x(2) - exp(-x(2)) ] ;
% 可以写成匿名函数的形式：
% h = @(x) [2 -1; -1 2]x - [exp(-x(1)); exp(-x(2))]
x0 = [ 3 3 ]; % test for x0 = [ 0 0 ] & [ -100 100 ]
options = optimset(‘Display’, ‘iter’);
% [ x, value ] = fsolve( @funfsolve, x0, options )
[ x, value, exitflag, output, jacobian ] = fsolve( @funfsolve, x0, options )
%% 例：求解满足XXX = [ 1 2; 3 4 ] 的矩阵X
h = @(x) xxx - [ 1 2; 3 4];
[x, value, flag ] = fsolve( h, [1 1; 1 1] )
%% 例：求解方程 x3 – 2x2 – x + 2 = 0
h = @(x) x.^3 - 2x.^2 - x + 2;
x = -2:0.1:3;
plot(x, h(x), ‘-r’);
grid on;
% 求解方法对比：
r1 = fsolve( h, 0 )
r2 = fzero ( h, 0 )
r3 = roots( [1 -2 -1 2 ] )
优化函数fzero/fminbnd/fminsearch/fsolve…小结：

在待求解域内，函数是连续的。求解时总是假设问题有解，但有些问题不收敛或收敛速度极慢，计算会中断。
只能得到一个解。如问题有多解，只能尝试改变不同的设定初值。
恰当设定初值，否则问题可能无解。
对多项式、线性方程组等问题不要用此类求解器。
更多的优化函数，参阅Optimization Toolbox。
直接使用交互式优化工具optimtool。

数据统计
常用函数
sum：求和
mean：平均值（数学期望）
median：中数
mode：出现最多的数值
max/min：最大值/最小值
var：方差（无偏估计量）(∑▒〖(x_i-x ̅)〗^2 )/(n-1)
std：标准差（无偏估计量）
cov：协方差（矩阵） (∑▒〖(x_i-x ̅ )(y_i-¯y)〗)/(n-1)
corrcoef：协方差系数（矩阵）

% 对多维数组，按列优先原则处理。
% 方差函数中，输入参数1为计算方差，无偏估计缺省为0。
% 协助方差函数只能按列处理，列数就是随机变量数目。
%% 例：单变量的方差和标准差
a = 1:10;
var(a) % 缺省：无偏估计
var(a,0) % 无偏估计/n-1
var(a,1) % 方差/n
v0 = sum( (a-mean(a)).^2 )/(length(a)-1)
v1 = sum( (a-mean(a)).^2 )/length(a)
std(a), std(a, 0), std(a, 1)
std(a)^2, std(a,1)^2 % 标准差的平方=方差
%% 多变量的计算
a = magic(3)
var(a)
cov(a) % 协方差矩阵
corrcoef(a)
a1 = a(:, 1); a2 = a(:, 2);
cov12 = sum( (a1-mean(a1)).*(a2-mean(a2)) ) / (length(a1)-1)
% 协方差如需按行计算，先做转置运算
cov(a’)
% 例：三个城市某月的日最高气温数据
t= [ 12 8 18
15 9 22
12 5 19
14 8 23
12 6 22
11 9 19
15 9 15
8 10 20
19 7 18
12 7 18
14 10 19
11 8 17
9 7 23
8 8 19
15 8 18
8 9 20
10 7 17
12 7 22
9 8 19
12 8 21
12 8 20
10 9 17
13 12 18
9 10 20
10 6 22
14 7 21
12 5 22
13 7 18
15 10 23
13 11 24
12 12 22 ];
scales = size( t )
plot( t ), legend(‘City-A’, ‘City-B’, ‘City-C’); xlabel(‘Days’); ylabel(‘Temperature’);

avg = mean(t) % 平均气温
dev = std(t) % 气温波动

% 每日气温与月平均值的差异
ft = t - avg(ones(length(t),1),

MATLAB-数值方法相关推荐

matlab期权定价模型比较,期权定价模型与数值方法（Matlab+Jupyter Notebook）
[实例简介] 期权定价模型与数值方法(Matlab+Jupyter Notebook) 这是在JupyterNotebook上运行的Matlab代码,其中包括:隐含波动率计算.二叉树模型.欧式期权蒙特 ...
二阶偏微分方程组龙格库塔法_数值方法（MATLAB版）（原书第3版）[Numerical Methods Using MATLAB，Third Edition]pdf...
摘要本书特点强大的图形表达宽泛的计算方法重点科学领域的重要算法大量可运行的实例数值方法(MATLAB版)(原书第3版)[Numerical Methods Using MATLAB,Thi ...
MATLAB秦九韶多项式求值算法的原理和迭代法求解的近似数值方法。
1..熟悉常用的Matlab操作: 2.了解秦九韶多项式求值算法的原理和迭代法求解的近似数值方法. 秦九韶多项式求值算法: 迭代法求解的近似数值: x=2; for k=1:10x=(x+2/x)/2 ...
数值方法的圣经-《应用数值方法（MATLAB实现）》第二版
Applied Numerical Methods Using MATLAB, 2nd Edition Applied Numerical Methods Using MATLAB 作者:Wenwu ...
牛顿科特斯型matlab,工程与科学数值方法的MATLAB实现(第4版)[PDF][119.29MB]
内容简介全书共分6大部分.第1部分介绍数值方法的背景知识.MATLAB的软件环境和编程模式,后5部分集中介绍数值方法的主要应用领域,具体包括求根与*化.线性代数方程组的求解.曲线拟合.数值积分与微分 ...
分数阶偏微分差分方程MATLAB,分数阶偏微分方程及其数值方法.ppt
分数阶偏微分方程及其数值方法郭柏灵北京应用物理与计算数学研究所定理6 设 (47) 为阶分数阶ODE,令 (48) 为相应的指数多项式,令 (49) 则如为最小整数, 为(47)的个线性无 ...
matlab数值解方程,[原创]MATLAB的Mupad应用之以数值方法解方程
本帖最后由 lijinfeng042 于 2010-11-3 19:26 编辑使用Mupad可以非常方便的求解非线性方程的数值解,简化输入的步骤,同时添加各种限制~只要有一下几个函数numeric: ...
用matlab画曲顶柱体费用数据,数值积分的matlab实现
实验10 数值积分实验目的: 1．了解数值积分的基本原理: 2．熟练掌握数值积分的MATLAB 实现: 3．会用数值积分方法解决一些实际问题. 实验内容: 积分是数学中的一个基本概念,在实际问题中也 ...
matlab平滑曲线_说说地震波的那些事儿（二）——地震影响曲线
原创:Eva哈哈哈转载自微信公众号:非解构做结构设计的朋友,一定很熟悉下面这张图.没错,它是抗规5.1.5条中的地震影响系数曲线.我们根据该曲线确定结构所受地震力的大小.作为结构抗震设计重要的一部 ...
matlab中表示拉普拉斯分布_CHAPT1：场论；电磁学和微波学的基本的数学手段和表示...
物理学中把某个物理量在空间一个区域内的分布称为场.从各种场的取值性质来看可以分成两大类,一类是每个点对应一个数值,这种场统称为标量场,如温度场.密度场等;另一类是每个点对应一个向量,这种场称为向量场 ...

MATLAB-数值方法

MATLAB-数值方法相关推荐

最新文章

热门文章