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), 😃 %…
figure, plot(ft), legend(‘City-A’, ‘City-B’, ‘City-C’);

% 出现频率最高的气温值
[n, f, c] = mode( t ) % n：频次最高的（最小）数值，f：最高频次
% c：频次最高的（全部）数值

% 各城市温度的相关性
corrcoef(t)
% 连续两天的温度变化量
dt = diff( t ) ;
plot( dt ), legend(‘City-A’, ‘City-B’, ‘City-C’);
mean(dt), std(dt), corrcoef(dt)
% 残缺数表的处理
% 读取数据，如data = xlsread(’.xlsx’)
% 生成残缺数表
t( 5, 1 ) = NaN; t( 10, 2 )= NaN; t(15, 2 ) = NaN;
plot( t );
mean(t), std(t) % nan运算结果不确定
%% 编程处理
% 对sum,mean等函数可以按列分别处理
t2 = t( :, 2); % 提取一列
t2( find( isnan(t2) ) ) = [ ]; % 剔除nan
size(t2)
sum(t2), var(t2)
%% Matlab提供的nan-簇函数，可以直接处理含缺数据
nansum( t ), nanvar( t )
% nanvar处理剔除nan后的一列数据。
% nancov剔除含nan的整行数据，保持矩阵结构。
% 验证nanvar和nancov的处理方法：
a = magic(4); b = a;
a(1,1) = NaN; c = a( 2:4, : );
nanvar(a), var(b), var©
nancov(a), cov(b), cov©
数据插值
三次样条插值方法：
给定数据系列(xi, yi), i = 1, 2, …, N，用(N-1)个三阶多项式，逐段逼近。
插值函数f的约束条件：
yi = f(xi)
f(xi+) = f(xi-)
f(xi+)’ = f(xi-)’
f(xi+)’’ = f(xi-)’’
边界区域f ‘’'相等
三次样条插值函数—spline
格式：
yy = spline( x, y, xx)
pp = spline( x, y )
% (x, y)是离散的观测数据（x即间断点），xx待插值点
% yy为插值，pp为插值的三次多项式
% y可以是矩阵或多维数组
x = linspace( 0, 2pi, 6), y = sin(x)
xx = linspace( 0, 2*pi ); % 插值结果与插值点数目无关
yy = spline( x, y, xx );
pp = spline(x, y)
% 用pp计算
y1 = ppval( pp, xx);
isequal( y1, yy)
% 图形显示插值的效果
plot( x, y, ‘o’, xx, yy, ‘-r’ )

% 参考sin(x)的图形
hold on
plot( xx, sin(xx), ‘-b’ )
legend( ‘data’, ‘spline’, ‘sin(x)’)
hold off
分段三次Hermite插值—pchip
% 例：
x = -3:3; y = [-1 -1 -1 0 1 1 1]; plot( x, y, ‘o’ )
xx = -3:0.01:3; yy = pchip(x ,y, xx);
hold on, plot( xx, yy, ‘-r’ ), hold off
% 比较spline插值
ys = spline(x, y, xx);
hold on
plot( xx, ys, ‘-k’ )
legend( ‘data’, ‘pchip’, ‘spline’, 2 )
hold off
%查看pchip和spline的特点和差异。
doc pchip
一维插值— interp1
% 比较不同插值方法
x = linspace(0,2pi, 6); y=sin(x); % (x, y)为离散数据
xx = linspace( 0, 2pi, 21 ); % 插值点数目
y1 = interp1( x, y, xx, ‘linear’ ); % ~ default
y2 = interp1( x, y, xx, ‘nearest’ );
y3 = interp1( x, y, xx, ‘cubic’ ); % ~ pchip
y4 = interp1( x, y, xx, ‘spline’ );
x0 = 0:0.1:2*pi; y0 = sin(x0);
plot( xx, y1, xx, y2, xx, y3, xx, y4, x0, y0);
legend( ‘linear’, ‘nearest’, ‘cubic’, ‘spline’, ‘sin(x)’ )
例：分析听力敏感度s和频率h关系
h = [ 20:10:100 200💯1000 1500 2000:1000:10000 ];
s = [ 76 66 59 54 49 46 43 40 38 22 14 9 6 3.5 2.5 1.4 0.7 0 -1 -3 -8 -7 -2 2 7 9 11 12 ];
semilogx(h, s, ‘-o’ ), grid on
% 用不同的插值方法估算h = 2500的s值
interp1( h, s, 2500 )
interp1( h, s, 2500, ‘linear’)
interp1( h, s, 2500, ‘nearest’)
interp1( h, s, 2500, ‘cubic’)
interp1( h, s, 2500, ‘spline’)
% 估算最小值（听力门限）
h1 = linspace(1000, 6000);
s1 = interp1(h, s, h1, ‘spline’); % spline vs. cubic
i = find( h>=1000 & h<=6000);
semilogx( h(i), s(i), ‘or’, h1, s1, ‘-b’ );
grid on
[ smin, i ] = min(s1)
hmin = h1( i )
二维插值—interp2
% 例：分析海洋深度的观测数据
x = 0:0.5:4;
y = 0:0.5:6;
z = [ 100.0 99.0 100.0 99.0 100.0 99.0 99.0 99.0 100.0
100.0 99.0 99.0 99.0 100.0 99.0 100.0 99.0 99.0
99.0 99.0 98.0 98.0 100.0 99.0 100.0 100.0 100.0
100.0 98.0 97.0 97.0 99.0 100.0 100.0 100.0 99.0
101.0 100.0 98.0 98.0 100.0 102.0 103.0 100.0 100.0
102.0 103.0 101.0 100.0 102.0 106.0 104.0 101.0 100.0
99.0 102.0 100.0 100.0 103.0 108.0 106.0 101.0 99.0
97.0 99.0 100.0 100.0 102.0 105.0 103.0 101.0 100.0
100.0 102.0 103.0 101.0 102.0 103.0 102.0 100.0 99.0
100.0 102.0 103.0 102.0 101.0 101.0 100.0 99.0 99.0
100.0 100.0 101.0 101.0 100.0 100.0 100.0 100.0 99.0
100.0 100.0 100.0 100.0 100.0 99.0 99.0 99.0 99.0
100.0 100.0 100.0 99.0 99.0 100.0 99.0 100.0 99.0 ];
% 绘制图形
size(x), size(y), size(z)
surf( x, y, z )
% 比较不同的插值方法估算（2.2, 3.3）处的深度。
interp2( x, y, z, 2.2, 3.3 )
interp2( x, y, z, 2.2, 3.3, ‘linear’ )
interp2( x, y, z, 2.2, 3.3, ‘nearest’ )
interp2( x, y, z, 2.2, 3.3, ‘cubic’ )
interp2( x, y, z, 2.2, 3.3, ‘spline’ )
% 网格的细化
[ xx, yy ] = meshgrid( 0:0.1:4, 0:0.1:6 );
zz = interp2( x, y, z, xx, yy, ‘spline’ );
surf( xx, yy, zz )
% 估计深度的峰值
zmax = max(max(zz))
[i, j] = find( zmax == zz )
xymax = [xx(i,j), yy(i,j)]
多维分散数据的插值—scatteredInterpolant
% 基于Delaunay triangulation处理离散数据点。
% 参阅delaunay和相关函数说明。
% Matlab早期定义的griddata和TriscatterInterp功能相似。
% --------------------------------------------
%% 构建2D离散的数据点
rng(0); % 设置随机数
x = rand(100, 1)*4 - 2;
y = rand(100, 1)*4 - 2;
z = x.*exp(-x.^2-y.2);
plot3(x, y, z, ‘o’)
view(2)
%% 定义插值运算器
F =scatteredInterpolant( x, y, z );
%% 输出网格和绘图
[xx, yy] = meshgrid(-2:0.2:2);
zz = F(xx, yy);
mesh(xx, yy, zz)
hold on, plot3(x, y, z, ‘o’), hold off
%% 编写代码，查看插值的效果
数据拟合
线性拟合：mldivide
非线性拟合：lsqcurvefit
数据预处理：smooth
Matlab拟合工具：Basic Fitting & cftool
最小二乘法（Least square method）：
对数据系列(xi, yi)，按函数y = f(x, a)进行拟合，要求残差平方和最小化。
ϵ^{2=∑▒〖[y_i-f(x_i,a)]〗}2
多元函数取得最小值的条件是：
〖∂ϵ〗^2/(∂a_i )=0
求解上述方程组，得到待估参数a
% 线性拟合
% 利用左除运算，求解线性方程组
%% 由x, y的测试数据
x = [0 0.3 0.8 1.1 1.6 2.3];
y = [0.6 0.67 1.01 1.35 1.47 1.25];
% 用函数 y= a_1+ a_2 e^x+ a_3 e^(-x) 进行数据拟合
x = x.’;
y = y.’;
plot( x, y, ‘o’ );
% 定义和求解线性方程组：M A = y，A = [a1; a2; a3]
M = [ ones(size(x)), exp(x), exp(-x) ]; % 系数矩阵
% M行数=数据点数目，列数=参数数目
A = M\y
% 超定线性方程组的最小二乘解
xx = 0:0.1:2.5;
yy = A(1) + A(2)*exp(xx) + A(3)*exp(-xx);
plot(x, y, ‘o’, xx, yy, ‘.-r’)
legend(‘original’, ‘fitted’)
%% 用二阶多项式拟合上述问题，并比较polyfit拟合结果
M = [ x.*x, x, ones(size(x)) ];
A = M\y;
xx = 0:0.1:2.5; yy = A(1)*xx.^2 + A(2)*xx + A(3);
plot(x, y, ‘o’, xx, yy, ‘.-r’)
p = polyfit(x, y, 2); % 直接的多项式拟合
A – p.’
% 非线性拟合：lsqcurvefit函数
% 调用格式
% [ a, norm2, residual, exitflag, … ] = lsqcurvefit( fun, a0, x, y, … )
% 按 y=asin(bx+c)e^(-dx) 进行参数拟合
x = 0:10;
y = [1.65 0.915 -0.72 -1.29 -0.417 0.259 0.446 0.361 -0.038 -0.39 -0.236 ];
plot(x, y, ‘–or’);
h = @(a,x) a(1).*sin(a(2).*x+a(3)).*exp(-a(4)*x); % 定义拟合函数
a0 = [1; 1; 1; 1]; % 迭代初值
options = optimset( ‘Display’, ‘iter’, ‘TolFun’, 1e-6); % 设置求解器选项
[a, norm2, residual, exitflag] = lsqcurvefit( h, a0, x, y, 0, 10, options )

xx = 0:0.1:10;
yy = h(a, xx);
hold on; plot(xx, yy); legend(‘original’, ‘fitted’); hold off
% 对下列数据按〖 y=a_1 e〗^(a_2 x) 进行参数拟合
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 ];
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 ];

plot(xdata, ydata, ‘or’);

h = @ (a, x) a(1)*exp(a(2)*x);
a0 = [1; -1];
a = lsqcurvefit( h, a0, xdata, ydata )
xx = 0:40; yy = h(a, xx);
hold on; plot(xx, yy); legend(‘original’, ‘fitted’); hold off;
% 按 ln(y)=〖ln(a〗_1)+a_2 x 对上述问题进行参数拟合
x = xdata.’; y = log(ydata.’);
M = [ones(size(x)), x ];
aa = M\y;
a2 = [exp(aa(1)); aa(2) ]
% 光滑化处理
% smooth函数调用格式：
% yy = smooth( x, y, span, method, …)
%% 例：
randn( ‘state’, 0 ); % ……
x = 0:0.005:1;
y = x.^2 + x + randn(size(x))/10;
plot( x, y, ‘-r’ );
% 去噪音
yy = smooth( x, y, 0.2, ‘loess’ ); % try other methods.
hold on; plot( x, yy.’, ‘LineWid’, 3); hold off;
% 数据拟合
p1 = polyfit( x, y, 2) %
p2 = polyfit( x, yy.’, 2) %
% 对比残差平方和
R1 = norm(y - polyval(p1, x)) % sum( (y - polyval(p1, x)).^2 ).^0.5
R2 = norm(yy.’ - polyval(p2, x))
Matlab拟合工具：

简单拟合：
图形窗口ToolsBasic Fitting
交互式曲线拟合：
cftool
%% 例1：用cftool分析全球平均气温变化
year = 1880:10:1980;
temperature = [ 0.01 0.02 0.03 0.04 0.06…
0.08 0.10 0.13 0.18 0.24 0.32 ];
% 分别用多项式和指数形式 y = a + b.exp( c.(x-d) )进行拟合。
% 1800~2100年间的气温变化。
%% 例2：分析美国人口的增长趋势
year = 1900:10:1990;
population = [75.995 91.972 105.711 123.203 131.669…
150.697 179.323 203.212 226.505 249.633];
数值积分
% 函数trapz按给定的间隔抽样、梯形分割近似计算。
% trapz可以直接处理离散数据。
% 函数trapz的计算方法
x = [ 0 1 3 4 ];
y = [ 0 1 1 2 ];
plot(x, y, ‘-or’)
trapz( x, y)
% 计算humps函数在[-1, 2]区间的积分值
x = linspace( -1, 2, 200); y = humps(x); plot( x, y, ‘.-’);
format long
x = linspace(-1, 2, 20); % 更多的抽样点n=20, 50…
y = humps(x);
area = trapz(x, y)
% 函数cumtrapz计算 ∫_(x_0)^x▒f(x)dx
x = linspace(-1, 2);
y = humps(x);
z = cumtrapz(x,y);
plotyy(x, y, x, z)
% 新版本中推荐使用integral和integral2代替quad*
% 函数quad和quadl采用任意间隔抽样和高阶估计。
format long
area = quad(@humps, -1, 2)
area = quadl(@humps, -1, 2)
% 例：空间曲线的长度
% x(t) = sin(2t), y(t) = cos(t), z = t; t = [0, 3pi]
t = 0:0.1:3pi;
plot3( sin(2t), cos(t), t )
h = @(t) sqrt( 4cos(2t).^2 + sin(t).^2 + 1 );
curvelength = quadl( h, 0, 3pi )
广义积分—quadgk
doc quadgk % 查看几种积分函数的比较说明
% 无界点
format long
quadgk( @(x) log(x), 0, 1)
quadgk( @(x) abs(x).^-0.5, -1, 1 )
% 无穷区间
quadgk( @(x) exp(-x), 0, inf )
quadgk( @(x) 1./(1+x.^2), -inf, inf )
% （特殊点）
quadgk( @(x) sin(x)./x, -1, 1 ) % 2*sinint(1)
quadgk( @(x) sin(1./x), 0, 1 ) % sin(1)-cosint(1)
多重积分—dblquad/triquad/quad2d
% 计算 函数sin(x)cos(y) + 1在[0, pi] & [-pi, pi] 区域的二重积分
h = @(x,y) sin(x).*cos(y) + 1;
[ x, y ] = meshgrid( 0:pi/20:pi, -pi:pi/20:pi );
z = h(x, y);
surf(x, y, z)

volume = dblquad(h, 0, pi, -pi, pi)

% 精确解
exact_solution = 2*pi^2;
err = abs(volume - exact_solution)/exact_solution
% 函数triplequad计算三重积分。
% 如二维积分区域不是矩形可以用函数quad2d

% 计算函数 [ (x+y)1/2(1+x+y)2 ]-1 在 [0, 1] & [0, 1-x] 的二重积分值

[x, y] = meshgrid( 0:0.1:1 );
z = 1./ ( sqrt(x+y) .* (1+x+y).^2 );
surf(x, y, z);

h = @(x, y) 1./ ( sqrt(x+y) .* (1+x+y).^2 );
ymax = @(x) 1-x;
volume = quad2d( h, 0, 1, 0, ymax )

% 精确解
exact_solution = pi/4 - 1/2
数值微分
% 一维差分diff
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 ];

diff (x)
diff (y)

y(2:end) - y(1:end-1)

dydx = diff(y)./diff(x)

plot(x, y, ‘.-’, x(1:end-1), dydx, ‘-or’ );
legend( ‘y(x)’, ‘dydx’ );
% 微分描述函数局部特性，对函数的波形变化非常敏感。
% 处理离散数据时，最好先对数据进行拟合，再进行微分。
p = polyfit(x, y, 2)
dp = polyder( p )
dydx2 = polyval(dp, x)
hold on
plot(x, dydx2, ‘–’)
legend( ‘y(x)’, ‘dydx’, ‘dydx2’ )
hold off
% 一阶微分的几种格式
向前差分：( yi+1 – yi ) / ( xi+1 – xi )
向后差分：( yi – yi-1 ) / ( xi – xi-1 )
中心差分：( yi+1 – yi-1) / ( xi+1 – xi-1 )
% 比较三种差分格式
step = 1; % try step = 0.5, 0.25, 0.1…
x = 0:step:5 ; y = sin(x) ;

x1 = x(1:end-1); dydx1 = diff( y ) / step ;
x2 = x(2:end) ; dydx2 = diff( y ) / step ;
x3 = x(2:end-1) ; dydx3 = ( y(3:end)-y(1:end-2) ) / 2 / step ;

x0 = 0:0.01:5; y0 = cos(x0);

plot( x0, y0, x1, dydx1, ‘.-r’, x2, dydx2, ‘.-g’, x3, dydx3, ‘-ok’ );
legend(‘精确值’, ‘向前差分’, ‘向后差分’, ‘中心差分’)
%% 梯度函数gradient
∇F(x,y,…)=∂F/∂x n_x+∂F/∂y n_y+⋯
% 在绘图中的应用：流线图
[x, y, z] = peaks(20);
dx = x(1,2) - x(1,1); dy = y(2, 1) - y(1, 1);
[dzdx, dzdy] = gradient( z, dx, dy );
c=contour(x ,y, z, 10); clabel©;
hold on
quiver(x, y, dzdx, dzdy);
hold off
% 另一个例子
[x, y] = meshgrid(-2:0.2:2);
z = x .* exp(-x.^2 - y.^2);
[u, v] = gradient(z, 0.2, 0.2);
c = contour(x, y, z, 10); clabel©;
hold on, quiver(x, y, u, v); hold off
% %离散Laplace算子del2

∇^2 F(x,y,…)=(∂^2 F)/(∂x^2 )+(∂^2 F)/(∂y^2 )+⋯
% 调用格式：
% L = del2( F, hx, hy, …)
% L = 2F / 2n，n~维数
% 曲面的表面曲率
[x, y, z] = peaks;
dx = x(1, 2) - x(1, 1); dy = y(2, 1) - y(1, 1);
L = del2( z, dx, dy );
surf( x, y, z, abs(L) )
shading interp
常微分方程
微分方程分类：
常微分方程：
初值问题
非刚性问题—ode23/ode45/ode113
刚性问题—ode15s/ode23s/ode23t/odetb
微分代数方程—ode23t/ode15s
全隐式微分方程—ode15i
边值问题—bvp4c/bvp5c
延时微分方程（DDE）—dde23/ddesd
偏微分方程：
一维抛物/椭圆型方程—pdepe
交互式PDE工具—pdetool
Matlab求解常微分方程的格式要求：

标准格式： y’ = f( t, y)
线性隐式格式： M y’ = f( t, y)
全隐式格式： f( t, y, y’ ) = 0
高阶常微分方程必须降阶成一阶微分方程组。
doc ode23 % 查看ode-系列函数的用法和差异
例：求解常微分方程
y’ = -y + t + 1
y(0) = 1
（该问题的精确解为：e-t + t）
% ----------------------------------------------------------------------
% 求解常微分方程的基本格式
% 程序调试完成后，可以将脚本改为一个M函数。

f = @ (t,y) -y+t+1; % 形式较复杂的方程最好写成M函数或子函数。

tspan = [ 0, 1 ]; % 求解域
y0 = 1; % 初始条件
options = odeset(‘RelTol’, 1e-3, ‘MaxStep’, 0.1); % 设置求解选项……
[t, y] = ode45(f, tspan, y0, options ); % 使用格式……

yy = t + exp(-t); % 精确解
% plot(t, y, ‘or’, t, yy, ‘.-b’) % 对比计算结果

hold on
plot(x, abs(y-yy), ‘-or’) % 查看计算误差
hold off

% 改变求解器、求解选项等查看计算误差
例：求解van der Pol方程

IC: y(0) = 2, y’(0) = 0;

降阶处理：令y1 = y, y2 = y1’

IC: y1(0) = 2, y2(0) = 0
% -------------------------------------------------------
% van de Pol方程：非刚性/刚性问题求解
% mu = 1 非刚性问题
% mu = 1000 刚性问题

mu = 1;

f = @(t, y) [ y(2); mu.*(1-y(1).^2).*y(2) - y(1) ];

tspan = [ 0, mu*20 ];
y0 = [ 2; 0 ];
[t, y] = ode15s( f, tspan, y0 );

plot( t, y(:,1), ‘.-r’ )

% 查看mu = 1~1000，比较ode45和ode15s的求解。
例：求解代数微分方程

IC: y (0) = [ 0.8; 0.1; * ]

% ------------------------------------------------------------------
% 求解DAE问题
h = @(t, y) [ -0.2*y(1) + y(2).y(3) + 0.3y(1).y(2)
2y(1).y(2) - 5y(2).y(3) - 2y(2).^2
y(1) + y(2) + y(3) - 1 ] ;

y0 = [0.8; 0.1; 0.1]; % 初值y3可以是估计值
tspan = [0, 10];

M = [ 1 0 0
0 1 0
0 0 0 ]; % 定义Mass矩阵
options = odeset( ‘Mass’, M, ‘RelTol’,1e-4 );

[t, y] = ode15s( h, tspan, y0, options);

plot(t, y, ‘.-’), legend(‘y1’, ‘y2’, ‘y3’, 3)
%% 将原方程y3消去，用ode45验证解法。
例：求解边值问题
Equ: y’’ - 3y’ + 2y = 0
BC: y(0) = 0, y(1) = e2 - e
（精确解：y = e2x – ex）

首先将方程降为一阶，令y1 = y, y2 = y1’
y1’ = y2
y2’ = 3y2 - 2y1
BC: y1(0) = 0, y1(1) = e2 - e
% ---------------------------------------------------------------------
% bvp4c求解BVP问题

hfun = @(x, y) [ y(2); 3y(2)-2y(1) ]; % 定义待求解方程
hbc = @(ya, yb) [ ya(1); yb(1)-exp(2)+exp(1) ]; % 定义边界条件
solinit = bvpinit( linspace(0,1,21), [0 1] ); % 初始化条件
sol = bvp4c( hfun, hbc, solinit ); % 输出格式是固定的

yy = exp(2*sol.x) - exp(sol.x);
plot( sol.x, yy, ‘.-b’, sol.x, sol.y(1,:), ‘or’ ) % 对比精确解
% err = sol.y(1,:) - yy; plot( sol.x, abs(err), ‘-or’) % 查看计算误差

%% 求解选项option = bvpset(……)

% 求解边值问题
Equ: y’’ + |y| = 0
BC: y(0) = 0, y(4) = -2
% 方程有两个解
% 不同的初始化条件的影响
% 求解Singular BVP问题
Equ: y’’ + 2/x * y’ + y^5 = 0
BC: y’(0) = 0, y(1) = sqrt(3)/2
% doc bvpset 查看使用格式
% 精确解：y = 1./sqrt( 1 + x.^2/3 );
例：一维稳态热扩散问题

对于导热系数为k厚度为L的一维墙体，x = 0处保持恒温T0，x = L处与环境对流换热（环境温度为Tf，对流换热系数为h）；墙体内含均匀分布的内热源q。该传热问题描述为：

BC:
取简化参数：L = 1, k = 1, h = 1, T0 = 0, Tf = 1, q = 2.
Equ: T’’ = -2
BC: T(0) = 0, T’(1) = 1 - T

% ----------------------------------------------------------------------------
% bvp4c求解一维稳态热扩散问题
% 第三类传热边界条件
hfun = @(x, T) [ T(2); -2 ]; % 降为一阶方程
hbc = @(Ta, Tb) [ Ta(1); Tb(1)+Tb(2)-1 ]; % 直接写入第三类边界条件
solinit = bvpinit( linspace(0,1,11), [0 0] ); % 初始化条件

sol = bvp4c( hfun, hbc, solinit );

TT = -sol.x.^2 + 2*sol.x; % 精确解 T = -x2 + 2x
plot( sol.x, TT, ‘.-b’, sol.x, sol.y(1,:), ‘or’ )
偏微分方程
一维椭圆/抛物型偏微分方程—pdepe

方程一般形式：

IC:
BC:
（下标x, t表示偏导数）
% 例：一维瞬态热扩散问题
Equ:
IC:
BC:
% thermaldiffusion.m
% -----------------------------------------------------------------------------------------------------
function thermaldiffusion
m = 0; % 0, 1, 2分别为直角坐标系、柱坐标、球坐标
x = linspace(0,1,11); % 空间网格数由用户确定
t = linspace(0,2,5); % 时间离散按ode15s，这里仅为输出结果的节点数

options = odeset( ‘RelTol’, 1e-3 ); % 精度主要取决于空间网格数
T = pdepe( m, @equ, @ic, @bc, x, t, options );

% surf(x,t,T); xlabel( ‘Distance’ ); ylabel( ‘Time’ ); zlabel( ‘Temperature’ );
figure; Texact = exp(-t)'sin(pix); surf(x, t, T-Texact); view(0, 0); %比较精确解

% 注意方程、初始条件、边界条件的定义格式
function [c, f, s] = equ(x, t, u, DuDx)
c = pi^2; f = DuDx; s = 0;

function u0 = ic(x)
u0 = sin(pi*x);

function [pl,ql,pr,qr] = bc(xl, ul, xr, ur, t)
pl = ul; ql = 0; pr = pi * exp(-t); qr = 1;
PDE Toolbox简介

交互式工具pdetool（求解二维问题）

命令窗口中键入pdetool
打开文件（求解新问题先选择Heat Transfer模型或……）
定义几何尺度（画图）
PDE选项：选择方程类型、参数设定
定义边界条件：参数设定
定义网格（画网格）
定义求解选项（计算精度、瞬态问题的计算时间）
计算……
计算结果绘图
计算结果较粗糙时，细化网格重新计算
保存计算结果（FEM网格信息，温度，动画）至工作区
保存PDE问题到M文件后退出

一些后处理：
将暂存在工作区的计算结果保持至文件
implay(M)播放动画
查看特定位置的温度、热流分布等

例：用pdetool求解瞬态二维热扩散问题

图示矩形区域长0.2m，宽0.1m；圆1区域含均匀加热源，圆心位于（0.05, 0.05），半径0.025；圆2区域为空区域，圆心位于（0.15, 0.05），半径0.025。矩形区域外边界绝热，圆2边界处保持等温T0 = 0，初始温度一致为T0，求解整个区域温度分布随时间的变化。

矩形域：k1 = 10 W/mK，ρ1 = 3000 kg/m3，c1 = 2000 J/kgK
圆1域：k2 = 100 W/mK，ρ1 = 5000 kg/m3，c2 = 3000 J/kgK，Q = 10000 W/m3

控制方程：
矩形区域：k_1 ∇^2 T_1=ρ_1 c_1 (∂T_1)/∂t
圆1区域：k_2 ∇^2 T_2+Q=ρ_2 c_2 (∂T_2)/∂t
边界条件：
矩形外边界：(∂T_1)/∂n=0
圆2边界：T_2=0
矩形-圆1边界：温度和热流连续
初始条件：
T(t=0)=T_0