现代科学运算-matlab语言与应用
东北大学 http://www.icourses.cn/home/ (薛定宇)
《高等应用数学问题的MATLAB求解(第四版)》
代码可在matlab r2016b 运行。
08 数据插值、函数逼近问题的计算机求解
08.01 一维插值问题
一维插值 interp1() 函数的调用格式为
插值方法:
Linear – 默认方法
Pchip – 三次函数插值
Nearest – 最近点插值
Spline – 样条插值,推荐使用
运用外推法,再区间[
] 外的点的插值
例8-1 数据插值
已知的数据点来自函数
生成的数据进行插值处理
得出较平滑的插值曲线
生成的数据显示
x = 0: .12: 1;
y = (x.^2 - 3*x + 5).*exp(-5*x).*sin(x);
plot(x, y, x, y, 'o');
不同的插值算法
x1 = 0: .02: 1;
y1 = interp1(x, y, x1);
y0 = (x1.^2 - 3*x1 + 5).*exp(-5*x1).*sin(x1);
y2 = interp1(x, y, x1, 'pchip');
y3 = interp1(x, y, x1, 'spline');
y4 = interp1(x, y, x1, 'nearest');
plot(x1, [y1', y2', y3', y4'], ':', x, y, 'o', x1, y0)
% 验证结果
[max(abs(y0(1:49) - y2(1:49))), max(abs(y0 - y3)), max(abs(y0 - y4))]
例8-2 徒手绘制规划曲线
用户可以在 x-y 平面内选择一些点
x轴排序,密集点插值,绘图
function sketcher(vis)
x = [];
y = [];
i = 1;
h = [];
axes(gca);
while 1
[x0,y0,but] = ginput(1);
if but == 1
x = [x, x0];
y = [y, y0];
h(i) = line(x0, y0);
set(h(i), 'Marker', 'o');
i = i + 1;
else
break
end
end
if nargin == 1
delete(h);
end
[x ii] = sort(x);
y = y(ii);
xx = [x(1):(x(end)-x(1))/100: x(end)];
yy = interp1(x, y, xx, 'spline');
line(xx, yy)
Lagrange插值算法及应用
生成一个插值x向量
插值算法为
matlab函数调用格式: y = lagrange(
)
例8-3 Runge现象
测试函数
生成数据,进行Lagrange插值
function y = lagrange(x0, y0, x)
ii = 1:length(x0);
y=zeros(size(x));
for i = ii
ij = find(ii ~= i);
y1 = 1;
for j = 1:length(ij)
y1 = y1.*(x-x0(ij(j)));
end
y = y+y1*y0(i)/prod(x0(i)-x0(ij));
end
x0 = -1+2*[0:10]/10;
y0 = 1./(1+25*x0.^2);
x = -1: .01:1;
y = lagrange(x0, y0, x);
ya = 1./ (1+25*x.^2);
plot(x, ya, x, y, ':');
调用interp1()函数
y1 = interp1(x0, y0, x, 'pchip');
y2 = interp1(x0, y0, x, 'spline');
plot(x, ya, x, y1, ':', x, y2, '--');
例8-4 人口预测
某省人口数据在census.xls给出
可以根据给出的数据预报未来的人口
X = xlsread('E:\\Work\\201806\\math_4th\\census.xls', 'B5:C67');
t = X(:, 1);
p = X(:, 2);
t0 = t(1:5:end);
p0 = p(1:5:end);
t1 = 1949:2015;
y = interp1(t0, p0, t1, 'spline', 'extrap');
plot(t, p, t1, y, t0, p0, 'o');
已知样本点的定积分计算
已知函数散点
function y = quadspln(x0, y0, a, b)
f = @(x) interp1(x0, y0, x, 'spline');
y = integral(f, a, b);
函数调用格式: I = quadspln(x0, y0, a, b)
例8-5 利用插值求数值积分
利用样条插值三分求解
用30个采样点求解:
x0 = 0: pi/30: pi;
y0 = sin(x0);
I = quadspln(x0, y0, 0, pi)
08.02 二维与高维插值问题
二维网格数据的插值问题
样本点是网格形式给出的
二维插值的函数
其中,
选项类似于 interp1,建议使用 spline 选项
这里的数据一定要采用网格数据
网格数据应该由 meshgrid函数生成
例8-7 网格数据的插值举例
根据下述函数生成一些稀疏的网格数据
由下面的函数可以生成样本点网格数据
进行各种插值拟合,并比较拟合结果
绘制已知数据的网格图:
[x, y] = meshgrid(-3: .6 : 3, -2: .4 : 2);
z = (x.^2 - 2*x) .* exp(-x.^2-y.^2-x.*y);
surf(x, y, z);
axis([-3, 3, -2, 2, -0.7, 1.5]);
不同插值算法的比较
默认插值算法进行插值:
[x1, y1] = meshgrid(-3: .2: 3, -2: .2: 2);
z1 = interp2(x, y, z, x1, y1);
surf(x1, y1, z1);
axis([-3, 3, -2, 2, -0.7, 1.5]);
三次插值和样条插值,并分析误差
[x, y] = meshgrid(-3: .6 : 3, -2: .4 : 2);
z = (x.^2 - 2*x) .* exp(-x.^2-y.^2-x.*y);
[x1, y1] = meshgrid(-3: .2: 3, -2: .2: 2);
z1 = interp2(x, y, z, x1, y1, 'cubic');
z2 = interp2(x, y, z, x1, y1, 'spline');
surf(x1, y1, z1);
figure;
surf(x1, y1, z2);
z = (x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
figure;
surf(x1, y1, abs(z-z1));
figure;
surf(x1, y1, abs(z-z2));
二维一般分布数据的插值问题
实际应用中能测到的并不总是网格分布的数据,经常需要解决样本点为散点的插值
griddata()函数的调用格式:z = griddata(
, ‘v4’)
其中,
是已知的二维样本点向量
x,y 是期望的插值位置
z表示插值的结果,维数和x,y一致
‘v4’是matlab 4.0 版本中提供的插值算法
例8-8 随机分布样本点的插值
给定函数
在矩形区域
内随机生成一组样本点
用griddata()进行插值处理, 并误差分析
数据生成与显示
x = -3 + 6*rand(200, 1);
y = -2 + 4*rand(200, 1);
z = (x.^2 - 2*x).*exp(-x.^2-y.^2-x.*y);
plot(x, y, 'x');
figure;
plot3(x, y, z, 'x');
不同插值方法比较
使用’cubic’和’v4’算法:
[x1, y1] = meshgrid(-3: .2: 3, -2: .2: 2);
z1 = griddata(x, y, z, x1, y1, 'cubic');
surf(x1, y1, z1);
z2 = griddata(x, y, z, x1, y1, 'v4');
figure;
surf(x1, y1, z2);
% 比较
figure;
z0 = (x1.^2 - 2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
surf(x1, y1, abs(z0-z1));
figure;
surf(x1, y1, abs(z0-z2));
例8-9 剔除某些样本点
给定的样本点在x-y平面分布较均匀,现在人为剔除某些点,表明已知数据分布不均匀,这时再进行插值分析,观察插值效果。
剔除圆内样本点
x = -3 + 6*rand(200, 1);
y = -2 + 4*rand(200, 1);
ii = find((x+1).^2 + (y+0.5).^2 > 0.5^2);
x = x(ii);
y = y(ii);
z = z(ii);
plot(x, y, 'x');
t = [0: .1: 2*pi, 2*pi];
x0 = -1 + 0.5*cos(t);
y0 = -0.5 + 0.5*sin(t);
line(x0, y0);
重新插值
% 用新样本点拟合曲面
[x1, y1] = meshgrid(-3: .2: 3, -2: .2: 2);
z1 = griddata(x, y, z, z1, y1, 'v4');
surf(x1, y1, z1);
% 误差分析
z0 = (x1.^2 - 2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
surf(x1, y1, abs(z0-z1)),
axis([-3, 3, -2, 2, 0, 0.1]);
% 误差等高线图
figure;
contour(x1, y1, abs(z0-z1), 30);
hold on;
plot(x, y, 'x');
hold off;
line(x0, y0)
高维插值问题
三维的网格数据生成:
[x, y, z] = meshgrid(
);
高维网格数据的生成
插值函数
网格数据插值 interpn
散点数据插值 griddatan
数据格式转换
ndgrid与meshgrid格式不同
% 格式转换方法
function [x1, y1, z1, v1] = mesh2nd(x, y, z, v)
switch nargin
case 3, x1 = x.'; y1 = y.'; z1 = z.';
case 4, z1 = z;
for i = 1: size(x, 3), x1(:, :, i) = x(:, :, i).';
y1(:, :, i) = y(:, :, i).';
v1(:, :, i) = v(:, :, i).';
end;
otherwise error('Error in input arguments')
end;
例8-10 三元函数插值
假设已知某三元函数
可以通过该函数生成一些网格样本点,是根据样本点进行拟合,并给出拟合误差
[x, y, z] = meshgrid(0: 0.3: 2);
[x0, y0, z0] = meshgrid(0: 0.1: 2);
V = sqrt(x.^x + y.^((x+y)/2) + z.^((x+y+z)/3));
V0 = sqrt(x0.^x0 + y0.^((x0+y0)/2) + z0.^((x0+y0+z0)/3));
V1 = interp3(x, y, z, V, x0, y0, z0, 'spline');
vol_visual4d(x0, y0, z0, V1);
08.03 样条插值问题
三次样条函数及其matlab表示
样本点
其中
S(x) 为三次样条函数的三个条件:
每个子区间
上, S(x)位三次多项式
S(x) 在整个区间 [x_i, x_{i+1}] 上有连续1-2阶导数
例8.12 正弦函数插值
稀疏数据的三次样条插值结果
x0 = [0, 0.4, 1 2, pi];
y0 = sin(x0);
sp = csapi(x0, y0);
fnplt(sp, ':');
hold on;
ezplot('sin(t)', [0, pi]);
plot(x0, y0, 'o');
多个自变量网格数据插值
网格数据三次样条插值类:
S = csapi(
, z)
为自变量的网格数据
为网格数据的样本点
得出的S是三次样条函数对象
得出的样条函数类可用下面函数处理
fnplt(), fnval()
例8-14 二维数据插值
给定
用三次样条插值方法得出网格数据的样条插值拟合,并绘制曲面
x0 = -3: .6: 3;
y0 = -2: .4: 2;
[x, y] = ndgrid(x0, y0);
z = (x.^2 - 2*x) .* exp(-x.^2 - y.^2 -x.*y);
sp = csapi({x0, y0}, z);
fnplt(sp);
B样条函数及其matlab表示
建立B样条函数对象
S = spapi(k, z, y)
其中k为用户选定的B样条阶次,一般选择k= 4, 5
例8-15 插值类的比较
给定正弦函数
生成稀疏数据
进行5次B样条插值函数拟合
与三次样条函数拟合的结果比较
x0 = [0, 0.4, 1 2, pi];
y0 = sin(x0);
ezplot('sin(t)', [0, pi]);
hold on;
sp1 = csapi(x0, y0);
fnplt(sp1, '--');
sp2 = spapi(5, x0, y0);
fnplt(sp2, ':');
hold off;
另一个函数
原函数
插值效果比较
x = 0: .12: 1;
y = (x.^2-3*x+5).*exp(-5*x).*sin(x);
ezplot('(x^2-3*x)*exp(-5*x)*sin(x)', [0, 1]);
hold on;
sp1 = csapi(x, y);
fnplt(sp1, '--');
sp2 = spapi(5, x, y);
fnplt(sp2, ':');
hold off;
08.04 基于样条插值的数值微积分
基于样条插值的数值微分运算
求单变量k阶导数的函数调用格式
= fnder(S, k)
多变量函数的偏导数的函数调用格式
= fnder(S, [
])
S可用三次或B样条类
例8-16数值微分
给定
生成一些数据点
用三次分段多项式样条函数与B样条插值函数,求出该函数的导数
syms x;
f = (x^2-3*x+5)*exp(-5*x)*sin(x);
ezplot(diff(f), [0, 1]);
hold on;
x = 0: .12: 1;
y = (x.^2-3*x+5).*exp(-5*x).*sin(x);
% 数值微分
figure;
sp1 = csapi(x, y);
dsp1 = fnder(sp1, 1);
fnplt(dsp1, '--');
sp2 = spapi(5, x, y);
dsp2 = fnder(sp2, 1);
fnplt(dsp2, ':');
hold off;
% 高阶导数,如三阶导数
figure;
ezplot(diff(f, 3), [0, 1]);
hold on;
dsp1 = fnder(sp1, 3);
fnplt(dsp1, '--');
dsp2 = fnder(sp2, 3);
fnplt(dsp2, ':');
hold off;
例8-17 二元函数
给定
生成一些数据点
利用数值插值的方法拟合曲面
并与解析解法绘制出的曲面相比较
x0 = -3: .3 : 3;
y0 = -2: .2: 2;
[x, y] = ndgrid(x0, y0);
z = (x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
sp = spapi({5, 5}, {x0, y0}, z);
dspxy = fnder(sp, [1, 1]);
fnplt(dspxy);
hold on;
syms x y;
z = (x^2-2*x)*exp(-x^2-y^2-x*y);
ezsurf(diff(diff(z, x), y), [-3, 3], [-2, 2])
样条插值的数值积分计算
数值积分运算
积分函数:
= fnint(S)
积分函数的结果为样条对数
在区间[a, b]求取定积分
I = diff(fnval(fnint(S), [a, b]))
例8-18 正弦函数的数值积分
给定
生成一些稀疏的样本点
用样条积分的方程求出定积分及积分函数
x = [0, 0.4, 1 2, pi];
y = sin(x);
sp1 = csapi(x, y);
a = fnint(sp1, 1);
xx = fnval(a, [0, pi]);
I1 = xx(2) - xx(1);
sp2 = spapi(5, x, y);
b = fnint(sp2, 1);
xx = fnval(b, [0, pi]);
I2 = xx(2) - xx(1);
% 积分函数绘制
ezplot('-cos(t)+2', [0, pi]);
hold on;
fnplt(a, '--');
fnplt(b, ':');
hold off;
基于样条插值的数值积分运算
数值积分运算
积分函数
= fnint(S)
积分函数的结果为样条函数
在区间[a, b]求取定积分
I = diff(fnval(fnint(S), [a, b]))
例8-19 双重积分的积分曲面
积分公式
fnint()函数不能求解,可以考虑fnder()函数
微分阶次选为负值可以计算积分
x0 = -2: 0.1 : 2;
y0 = -1: 0.1: 1;
[x, y] = ndgrid(x0, y0);
z = exp(-x.^2/2).*sin(x.^2+y);
S = spapi({5, 5}, {x0, y0}, z);
S1 = fnder(S, [-1, -1]);
S2 = fnval(S1, {x0, y0});
surf(y0, x0, S2);
08.05 由已知数据拟合函数
多项式拟合
多项式拟合就是用多项式拟合已知数据
多项式拟合函数调用格式
p = polyfit(x, y, n)
其中:
x和y为原始的样本点构成的向量
n为选定的多项式阶次
p为多项式系数按降幂排列得出的行向量
例8-20 数据的多项式拟合
原型函数
使用该函数生成一些数据点
用多项式方法在不同的阶次下进行拟合
拟合该数据的3次多项式
x0 = 0: 0.1: 1;
y0 = (x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);
p3 = polyfit(x0, y0, 3);
x = 0: .01: 1;
ya = (x.^2-3*x+5).*exp(-5*x).*sin(x);
y1 = polyval(p3, x);
plot(x, y1, x, ya, x0, y0, 'o');
不同阶次拟合
p4 = polyfit(x0, y0, 4);
y2 = polyval(p4, x);
p5 = polyfit(x0, y0, 5);
y3 = polyval(p5, x);
p8 = polyfit(x0, y0, 8);
y4 = polyval(p8, x);
plot(x, ya, x0, y0, 'o', x, y2, x, y3, x, y4);
% Taylor幂级数展开
syms x;
y = (x^2-3*x+5)*exp(-5*x)*sin(x);
p = taylor(y, 'Order', 9)
例8-21 数据拟合
对下式进行多项式拟合,并观察拟合效果
就不同的多项式阶次进行曲线拟合:
x0 = -1+2*[0:10]/10;
y0 = 1./(1+25*x0.^2);
x = -1: 0.01: 1;
ya = 1./(1+25*x.^2);
p3 = polyfit(x0, y0, 3);
y1 = polyval(p3, x);
p5 = polyfit(x0, y0, 5);
y2 = polyval(p5, x);
p8 = polyfit(x0, y0, 8);
y3 = polyval(p8, x);
p10 = polyfit(x0, y0, 10);
y4 = polyval(p10, x);
plot(x, ya, x, y1, x, y2, x, y3, x, y4);
函数线性组合的曲线拟合方法
已知函数的线性组合
Ac = y
例8-21 多项式拟合
原型函数
进行多项式拟合的各个函数为
观察多项式拟合效果
x = [0: .1: 2]';
n = 7;
A = [];
y = (x.^2+3*x+5).*exp(-5*x).*sin(x);
for i = 1:n+1
A(:, i) = x.^(n+1-i);
end;
c = A\y
08.06 最小二乘曲线拟合
最小二乘曲线拟合
给定一组数据
已知函数原型
最小二乘曲线拟合的目标为求出一组待定系数的值,使得目标函数为最小
matlab求解
[a,
] = lsqcurvefit(Fun,
, x, y, options)
例8-25 最小二乘曲线拟合
由下面的语句生成一组数据
x = 0: 0.1: 10;
y = 0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);
该数据满足原型函数:
% 最小二乘拟合
f = @(a, x)a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x);
[xx, res] = lsqcurvefit(f, [1, 1, 1, 1, 1], x, y)
例8-26 最小二乘拟合
已知数据可能满足
求满足数据的最小二乘解 a, b, c, d的值
令
原型函数
% 输入已知数据:
x = 0.1: 0.1: 1;
y = [2.3201, 2.6470, 2.9707, 3.2885, 3.6008, ...
3.9090, 4.2147, 4.5191, 4.8232, 5.1275];
% 求解
f = @(a, x)a(1)*x + a(2)*x.^2 .* exp(-a(3)*x) + a(4);
a = [1;2;2;3];
[a, b, c, d] = lsqcurvefit(f, a, x, y);
a, b
多变量函数的最小二乘函数拟合
多变量待定函数
目标函数:
最小二乘拟合
例8-27 多变量函数拟合
原型函数
数据文件:c8data1.dat
令
f = @(a, X)a(1)*X(:, 1) .^ (a(2)*X(:, 1)) + ...
a(3)*X(:, 2) .^ (a(4)*(X(:, 1) +X(:, 2))) + ...
a(5)*X(:, 3) .^ (a(6)*(X(:, 1) + X(:, 2) + X(:, 3)));
% XX = load('E:\\Work\\201806\\math_4th\\c8data1.dat');
XX = load('c8data1.dat');
X = XX(:, 1:3);
v = XX(:, 4);
a0 = [2 3 2 1 2 3];
[a, fm] = lsqcurvefit(f, a0, X, v)
08.07 函数的有理式逼近
连分式的数学形式
调用格式:cf = contfrac(f, n)
[cf, r] = contfrac(f, n, ‘x=a’)
function key=isanumber(a)
key=0; if length(a)~=1, return; end
switch class(a)
case 'double', key=1;
case 'sym', try, double(a); key=1; catch, end
end
function [cf, r] = contfrac(f, varargin)
[n, a] = default_vals({6,0}, varargin{:});
if isanumber(f)
cf = feval(symengine, 'contfrac', f, n);
p1 = char(cf);
k = strfind(p1, ',');
k1 = strfind(p1, '/');
if nargout > 1
r = sym(p1(k(end)+1:k1-1))/sym(p1(k1+1:end-1));
end
else
if isfinite(a)
str=num2str(a);
else
str='infinity';
end
cf = feval(symengine, 'contfrac', f, ['x=' str], n);
if nargout > 1
r = feval(symengine, 'contfrac::rational', cf);
end
end
例8-28
的连分数近似
连分数近似cf = contfrac(pi, 20)
结果: 292大于其它值
截断于该项103993/33102
近似值: 3.141592653012
例8-29 函数的连分式近似
原型函数
10级连分式
syms x;
f = sin(x)*exp(-x)/(x+1)^3;
[cf, y] = contfrac(f, 10)
% 不同连分式近似
[cf1, f8] = contfrac(f, 8),
[cf2, f10] = contfrac(f, 10)
连分式拟合效果
% 拟合效果
ezplot(f, [0, 2]);
hold on;
ezplot(f8, [0, 2]);
ezplot(f10, [0, 2]);
figure;
ezplot(f, [0, 5]);
hold on;
ezplot(f8, [0, 5]);
ezplot(f10, [0, 5]);
% 参考点x = 0.5
[cf1, f8] = contfrac(f, 8, 0.5)
ezplot(f, [0, 5]);
hold on;
ezplot(f8, [0, 5]);
hold off;
有理式拟合:Pad
近似
函数 f(s) 的幂级数展开
系数计算
令相应的 s 项系数的值相等,得出
其中:
且有
编写函数计算
给定函数的
有理数近似调用格式
[n, d] = padefcn(c, r, m)
function [nP, dP] = padefcn(c, r, m)
w = -c(r+2:m+r+1)';
vv = [c(r+1:-1:1)';
zeros(m-1-r,1)];
W = rot90(hankel(c(m+r:-1:r+1),vv));
V = rot90(hankel(c(r:-1:1)));
x = [1 (W\w)'];
y = [1 x(2:r+1)*V'+c(2:r+1)];
dP = x(m+1:-1:1)/x(m+1);
nP = y(r+1:-1:1)/x(m+1);
例8-30 函数的有理近似
给定函数
求有理函数近似
分别选择分子阶次为0和不同的分母阶次
syms x;
c = taylor(exp(-2*x), 'Order', 10);
c = sym2poly(c);
c = c(end:-1:1);
ezplot(exp(-2*x), [0, 4]);
x = 0: 0.01: 8;
nd = [3:7];
for i = 1: length(nd)
[n, d] = padefcn(c, 0, nd(i));
y = polyval(n, x)./polyval(d, x);
line(x, y);
end;
符号运算
符号运算版本的padefcnsym()函数
function c = polycoef(p, x)
c=[];
n = 0;
p1=p;
n1=1;
nn=1;
if nargin == 1
x=symvar(p);
end
while (1)
c = [c subs(p1, x, 0)];
p1 = diff(p1, x);
n = n+1;
n1 = n1*n;
nn = [nn, n1];
if p1 == 0
c=c./nn(1:end-1);
c=c(end:-1:1);
break;
end
end
function G = padefcnsym(f, r, m)
c = taylor(f, 'Order', r+m+1);
c = polycoef(c);
c = c(end: -1: 1);
w = -c(r+2:m+r+1)';
vv = [c(r+1:-1:1)';
zeros(m-1-r, 1)];
W = rot90(hankel(c(m+r:-1:r+1), vv));
V = rot90(hankel(c(r:-1:1)));
X = [1 (W\w)'];
y = [1 X(2:r+1)*V'+c(2:r+1)];
dP = X(m+1:-1:1)/X(m+1);
nP = y(r+1:-1:1)/X(m+1);
syms x;
G = poly2sym(nP, x)/poly2sym(dP, x);
例8-31 符号运算求取Pade近似
给定函数
用符号运算方式求取 Pade近似模型
syms x;
f = exp(-2*x);
G = padefcnsym(f, 0, 8)
最终结果
08.08 特殊函数
误差函数与补误差函数
误差函数
补误差函数
两个函数的关系: erfc(z) = 1 - erf(z)
例8-34 误差函数与补误差函数的曲线
直接绘制:y = erf(z), y = erfc(z)
x = -5: 0.1: 5;
y1 = erf(x);
y2 = erfc(x);
plot(x, y1, x, y2, '--');
Gamma函数
Gamma函数的数学描述
Gamma函数调用格式: y = gamma(x)
其中,x可以为向量,矩阵或其它数据结构
gamma()函数可用于数据与符号运算
特殊性质:
常用性质:
例8-37 Gamma函数曲线
试绘制(-4, 3)区间内Gamma函数的曲线
x = -4: 0.0001: 3;
plot(x, gamma(x));
ylim([-40, 40]);
注意:
例8-38 不可积函数的解析求解
试写出下面无穷积分的“解析解”
syms m n t z;
syms x positive;
I1 = int(sin(t)^(2*m-1)*cos(t)^(2*n-1), t, 0, pi/2)
I2 = simplify(int(t^(x-1)*cos(t), t, 0, inf))
不完整Gamma函数
不完整
-函数的数学描述
求取不完整
-函数的调用格式
y = gammainc(x, b)
Beta函数
Beta函数的数学描述
Beta函数的调用格式: y = beta(m, z)
Bessel函数
对于Bessel微分方程
该方程的通解为
Bessel函数计算
第一类Bessel函数调用格式: y = besselj(
)
第二类Bessel函数调用格式: y = bessely(
)
第三类Bessel函数调用格式: y = besselh(
)
其中,
为阶次
Legendre函数
对于Legendre微分方程
该方程的通解为
Legendre函数计算
计算关联Legendre函数
调用格式: Y = legendre(n, x)
其中, Y为一个矩阵,其各行分别为
该函数要求:
例8-41 Legendre 函数曲线
x = -1: 0.04: 1;
Y = legendre(2, x);
plot(x, Y);
超几何函数
超几何函数的一般形式
符号运算: y = hypergeom(
)
数值运算:
= subs(
)
08.09 Mittag-Leffler函数
Mittag-Leffler函数
指数函数的幂级数展开
单参数 Mittag-Leffler 函数
其它常用的单参数 Mittag-Leffler 函数
双参数 Mittag-Leffler 函数
双参数Mittag-Leffler函数的一些公式
% Mittag-Leffler解析函数
function f = mittag_leffler(aa, z)
aa = [aa 1];
a = aa(1);
b = aa(2);
syms k;
f = simplify(symsum(z^k/gamma(a*k+b), k, 0, inf));
函数入口:
f1 = mittag_leffler(
, z)
f1 = mittag_leffler([
], z)
例8-36 Mittag-Leffler函数
单参数 Mittag-Leffler 函数推导
syms z;
I1 = mittag_leffler(1/sym(3), z),
I2 = mittag_leffler(3, z)
I3 = mittag_leffler(4, z)
I4 = mittag_leffler(5, z)
% 在matlab2008a下可以求出非超几何函数
例8-37 双参数 Mittag-Leffler 函数推导
syms z;
I5 = mittag_leffler([4, 1], z),
I6 = mittag_leffler([4, 5], z),
I7 = mittag_leffler([5, 6], z),
I8 = mittag_leffler([1/sym(2), 4], z)
双参数 Mittag-Leffler 函数的整数阶导数
调用格式:
f = ml_func(
)
f = ml_func(
)
function f = ml_func(aa, z, varargin)
aa=[aa,1,1,1]; a=aa(1); b=aa(2); c=aa(3); q=aa(4); f=0; k=0; fa=1;
[n,eps0]=default_vals({0,eps},varargin{:});
if n==0
while norm(fa,1)>=eps0
fa=gamma(k*q+c)/gamma(c)/gamma(k+1)/gamma(a*k+b) *z.^k;
f=f+fa; k=k+1;
end
if ~isfinite(f(1))
if c*q==1
f=mlf(a,b,z,round(-log10(eps0))); f=reshape(f,size(z));
else, error('Error: truncation method failed'); end, end
else
aa(2)=aa(2)+n*aa(1); aa(3)=aa(3)+aa(4)*n;
f=gamma(q*n+c)/gamma(c)*ml_func(aa,z,0,eps0);
end
例8-38 Mittag-Leffler函数
得出
的解析解与数值解(比较耗时)
syms z;
I8 = mittag_leffler([1/sym(2), 4], z);
t = 0: 0.01: 2;
y = subs(I8, z, t);
y1 = ml_func([1/2, 4], t);
plot(t, y, t, y1, '--');
% 数值解
t = 0: 0.01: 0.5;
y1 = ml_func([1/2, 4], t);
plot(t, y1);
例8-46 绘制函数
假设 z 是复数变量
% 生成网格, 构造复数矩阵, 绘制实部表面图
[x, y] = meshgrid(-6: 0.2: 6);
z = x + sqrt(-1)*y;
L = ml_func([0.8, 0.9], z);
L1 = real(L);
ii = find(L1 > 3);
L1(ii) = 3;
ii = find(L1 < -3);
L1(ii) = -3;
surf(x, y, L1);
axis([-6 6 -6 6 -3 3])
% 虚部表面图
figure;
L2 = imag(L);
ii = find(L2 > 3);
L2(ii) = 3;
ii = find(L2 < -3);
L2(ii) = -3;
surf(x, y, L2);
axis([-6 6 -6 6 -3 3]);
多参数Mittag-Leffler函数
三参数Mittag_leffler函数
四参数Mittag_leffler函数
Pochhammer符号
Mittag-Leffler函数的导数
n阶导数公式
matlab求解:
f = ml_func(
)
f = ml_func(
)
f = ml_func(
)
f = ml_func(
)
08.10 信号的相关分析
信号的相关分析
信号
的自相关函数的定义:
且相关函数为偶函数,即
两个信号
的互相关函数为:
例8-39 自相关函数计算
给定信号
求它的自相关函数
syms A1 A2 w1 w2 t1 t2 t tau T;
x = A1*cos(w1*t + t1) + A2*cos(w2*t + t2);
Rxx = simplify(limit(int(x*subs(x, t, t+tau), t, -T, T)/2/T, T, inf))
结果:
相关分析的数值方法
在实验中测出两组数值:
计算两组数据的相关系数
求取已知向量的相关系数矩阵
R = corrcoef(x, y)
R = corrcoef([x, y])
例 8-40 相关系数矩阵
用下列两个原型函数分别生成一组数据, 并求取其相关系数矩阵
x = 0: 0.01: 5;
y1 = x.* exp(-4*x) .* sin(3*x);
y2 = x .* exp(-4*x) .* cos(3*x);
R = corrcoef(y1, y2)
相关函数计算
在实验中测出两组数据
定义 x 序列的自相关函数
其中, m < n, 自相关函数是偶函数
定义出相关函数
求自相关的调用格式
= xcorr(x, N)
求互相关的调用格式
= xcorr(x, y, N)
其中, N为k的最大取值, 可以忽略
例8.41 相关函数
用函数
生成一组数据,用数值方法计算相关函数
并和已知理论曲线进行比较
t = 0: 0.01: 5;
x = t.*exp(-4*t).*sin(3*t);
y = t.*exp(-4*t).*cos(3*t);
N = 150;
c1 = xcorr(x, N);
x1 = [-N: N];
stem(x1, c1);
figure;
c1 = xcorr(x, y, N);
x1 = [-N: N];
stem(x1, c1)
08.11 信号滤波与滤波器设计
例8-42 污染的信号
给定信号
叠加标准误差为
的零均值的白噪声信号
绘制噪声污染后的信号曲线
x = 0: 0.002: 2;
y = exp(-x).*sin(5*x);
r = 0.05*randn(size(x));
y1 = y + r;
plot(x, y1);
线性滤波器的一般模型
离散线性时不变滤波器模型一般表达式为:
假设输入信号为 x(n), 则经过滤波器后的输出信号为:
三种滤波器类型
根据n和m的不同取值,定义常用的滤波器
FIR滤波器:又称为有限长脉冲响应滤波器, m = 0, a为标量,称为移动平均模型(MA)
IIR滤波器:又称全极点无限长脉冲响应滤波器,也称自回归模型(AR),n = 0, b为标量
ARMA滤波器:又称为一般IIR滤波器,也称为自回归移动平均模型(ARMA),要求n与m均不为0
滤波的信号向量调用格式: y = filter(b, a, x)
a,b 向量表示滤波器, x为需要过滤的信号
对滤波器进行放大倍数分析
[h, w] = freqz(b, a, N)
只获得放大倍数的幅值
plot(w, abs(h))
例8-43 滤波器效果
假设已经设计滤波器H(z)如下, 试得出该滤波器的放大倍数, 并观察滤波后的信号
b = 1.2296e-6*conv([1 4 6 4 1], [1 3 3 1]);
a = conv([1, -0.7265], ...
conv([1, -1.488, 0.5644], ...
conv([1, -1.595, 0.6769], ...
[1, -1.78, 0.8713])));
x = 0: 0.002: 2;
y = exp(-x).*sin(5*x);
r = 0.05*randn(size(x));
y1 = y+r;
[h, w] = freqz(b, a, 100);
plot(w, abs(h))
figure;
y2 = filter(b, a, y1);
plot(x, y1, x, y2);
滤波器及matlab实现
常用的滤波器类型
Butterworth滤波器
[b, a] = butter(n,
)
Chebyshev I 型滤波器
[b, a] = cheby1(n, r,
)
Chebyshev II型滤波器
[b, a] = cheby2(n, r,
)
例8-44 滤波器仿真
给定信号
叠加标准差为
的零均值的白噪声信号
试对阶次及不同的
值, 设计Butterworth 滤波器并比较滤波效果
f1 = figure;
f2 = figure;
for n = 5: 2: 20
figure(f1);
[b, a] = butter(n, 0.1);
y2 = filter(b, a, y1);
plot(x, y2);
hold on;
figure(f2);
[h, w] = freqz(b, a, 100);
plot(w, abs(h));
hold on;
end;
阶次为7时
f1 = figure;
f2 = figure;
for wn = 0.1: 0.1: 0.7
figure(f1);
[b, a] = butter(7, wn);
y2 = filter(b, a, y1);
plot(x, y2);
hold on;
figure(f2);
[h, w] = freqz(b, a, 100);
plot(w, abs(h));
hold on;
end
其它滤波器
设计高通滤波器
[b, a] = butter(n,
, ‘high’)
设计带通滤波器
[b, a] = butter(n, [
])