matlab数值计算

最新推荐文章于 2022-08-03 16:33:51 发布

mmk27

最新推荐文章于 2022-08-03 16:33:51 发布

阅读量3.5k

点赞数

分类专栏：数学建模

本文链接：https://blog.csdn.net/mmk27_word/article/details/104065871

版权

数学建模专栏收录该内容

6 篇文章 0 订阅

订阅专栏

第一节多项式运算

一、多项式的表示、求值、求根

1. 表示： matlab中把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。

p=[2,3,-1,4,5]表示多项式：

显示多项式的数学形式： p1=poly2str(p,'x') 如：p=[1,2,3,4] p1=poly2str(p,'x') 结果：p1 = x^3 + 2 x^2 + 3 x + 4

2. 求多项式在某点的值的函数：polyval(p,x) p：多项式对应的行向量, x：待求值的点

在x=1,2,-3,5点的值

p=[2,3,-1,4,5] x=[1,2,-3,5] y=polyval(p,x)

3 . 求多项式p的的根 m=roots(p)

如：p=[-1,2,3,1] x=roots(p)

二、多项式的运算

1. 四则运算

加减法：p1+p2 ； p1-p2

乘法：用卷积函数conv来实现如：p1=[1,0,-2,-5]; p2=[2,3,1]; conv(p1,p2)

除法：用deconv来实现对于任意两个多项式p1, p2, deconv(p1,p2)的值为两个行向量, 即[q,r]=deconv(p1,p2), 其中q是p1除以p2的商, r是余, 它们满足p1=conv(p2,q)+r.

如：p1=[1,0,-2,-5]; p2=[2,1]; [q,r]=deconv(p1,p2)

2. 求导、积分运算

求多项式p的导数：polyder(p)

求多项式p的不定积分：polyint(p) polyint(p,a,b) 非定积分

polyder(a,b): 求多项式a,b乘积的导数

polyint(a,b): 求多项式a,b乘积的积分

3. 特征多项式运算

（1）当x为行向量时，poly(x) 构造以x为根的多项式如：x=[-1,2,3] q=poly(x)

（2）当A为方阵时， poly(A) 构造以A的特征多项式。

第二节多项式拟合

一、数据拟合：

问题的提法：已知曲线上n个点(xi,yi), (i=1,2,…,n), xi互不相同，寻求一个函数y=F(x), 使y=F(x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。

基本思路:对n个点构成的数据集｛(xi,yi), i=1,2,…,n｝, m为小于n的整数，设是m+1个线性无关的连续函数（称为基本拟合函数），取按准则确定常数（最小二乘法）

按此方法得到的函数称为数据集(xi,yi), i=1,2,…,n的最小二乘拟合函数。

衡量拟合情况优劣的指标有：

（1）SSE（误差平方和）The sum of squares due to error 该统计参数计算的是拟合数据和原始数据对应点的误差的平方和：

SSE越接近于0，说明模型选择和拟合效果好，数据预测也成功。下面的指标MSE和RMSE与指标SSE有关联，它们的校验效果是一样的。

（2）MSE（方差）Mean squared error

（3）RMSE（剩余标准差）Root mean squared error 该统计参数，也叫回归系统的拟合标准差，是MSE的平方根

（4）R-square（判断系数，拟合优度）Coefficient of determination，判断系数是由两个参数SSR和SST决定的，对总平方和进行分解有SST =SSE +SSR，其中SSE 是误差平方和，反映随机误差对 y 的影响， SSR 称为回归平方和，反映自变量对y 的影响。判断系数定义为若它接近1，表示F(x)能较好的拟合该数据集。

（5）调整的判断系数统计学家主张在回归建模时，应采用尽可能少的自变量，不要盲目地追求判定系数 R2 的提高。其实，当变量增加时，残差项的自由度就会减少。而自由度越小，数据的统计趋势就越不容易显现。为此，又定义一个调整判定系数：

当 n 很大、 m 很少时，两者之间的差别不是很大；但是，当 n 较少，而m 又较大时，就会远小于

二、多项式拟合

若基本拟合函数为：

即用m次多项式拟合给定数据，称多项式拟合。

Matlab中的多项式拟合函数：polyfit

格式：a=polyfit(x,y,n) x:横坐标 y:纵坐标

功能：对一组点(x,y)进行n次多项式拟合。

x,y为要拟合的数据， n为拟合多项式的次数，

a为拟合多项式系数

拟合多项式在x处的值可用polyval(a,x)计算

可进入拟合工具箱cftool(x,y)处理数据拟合问题。例1 对x=0:0.1:1; y=[-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5,35,9.22] 求三次拟合多项式。

解   x=0:0.1:1;
       y=[-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];
       p=polyfit(x,y,3)
       xx=0:0.01:1;
      yy=polyval(p,xx);
      plot(xx,yy,'-b',x,y,'.r','markersize',20)

t=0:1:15;  
y=[30.0 29.1 28.8 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8];
 plot(t,y,'r*')

A=polyfit(t,y,3) 
 tt=0:0.5:15;
z=polyval(A,tt); 
plot(t,y,'r*',tt,z,'b')

对数据求出合适次数的拟合多项式p(x)，由弧长公式：得对河底光缆长度的估计s，所需光缆为s+y(0)+y(20) 6次多项式拟合较为合适。光缆总长：46.28m

x=0:1:20;
y=[-9.2,-8.91,-8.02,-7.95,-8.11,-9.05,-10.12,-11.08,-12.25,-13.26,-13.45,-12.62,-11.28,-10.23,-9.25,-7.91,-8.05,-8.86,-9.85,-10.81,-10.96]; 
polytool(x,y,3)
n=input('拟合次数n=')
format long
p=polyfit(x,y,n)
dp=polyder(p)
x1=0:0.01:20;
y1=sqrt(1+polyval(dp,x1).^2);
s=trapz(x1,y1)-y(1)-y(end)

第三节非线性拟合

问题的提法：设非线性函数其中是未知参数由实验获得该函数数据点，确定参数。

方法：最小二乘法

1. 非线性拟合函数：lsqcurvefit

格式：[x, resnorm,r,flag]=lsqcurvefit(fun, x0,xdata,ydata) x0为初始解向量；xdata，ydata为满足关系ydata=F(x, xdata)的数据； fun为待拟合函数， resnorm=sum ((fun(x,xdata)-ydata).^2)，即在x处残差的平方和； r=fun(x,xdata)-ydata，即在x处的残差； flag为终止迭代的条件

例1 确定模型中的参数，已知数据点：

xdata = [3.6,7.7,9.3,4.1,8.6,2.8,1.3,7.9,10.0,5.4];

ydata = [16,150.5,260.1,22.5,206.5,9.9,2.7,165.5,325.0,54.5];

解法：定义函数 f=@(a,x) a(1)*x.^2 + a(2)*x.*sin(x) + a(3)*x.^3;

由 [a, resnorm,r]=lsqcurvefit(f, a0,xdata,ydata)确定参数，画图观察拟合效果

2. 用fittype函数做曲线拟合

使用fittype函数可以自定义拟合函数，可以满足线性拟合和非线性拟合。fittype函数具有很灵活的配置，基本满足各种复杂场景。使用方法：

（1）做多项式拟合

p=fittype('poly2')

f=fit(x,y,p)

plot(f,x,y)

如：对下列数据做多项式拟合 x=[1;1.5;2;2.5;3]; y=[0.9;1.7;2.2;2.6;3];

x=[1;1.5;2;2.5;3];
y=[0.9;1.7;2.2;2.6;3];
p=fittype('poly2')
f=fit(x,y,p)
plot(f,x,y)

（2）做非线性拟合

方法1：（线性化方式）

如： ex= {'x.^2', ‘sin(x)','1'}; f=fittype(ex)

实验：x=[0.2,0.5,0.8,1.1,1.2,1.5,1.8,2]; y=[2.35,1.38,0.81,0.62,0.78,1.43,2.25,3.18];

方法2：（非线性化方式）如：p=fittype(‘a*x.^2+b*sin(x)+c’,‘independent’,‘x’)

由物理背景知：

 clc,clear
x=[0,0.4,1.2,2,2.8,3.6,4.4,5.2,6,7.2,8,9.2,10.4,11.6,12.4,13.6,14.4,15]';
y=[1,0.85,0.29,-0.27,-0.53,-0.4,-0.12,0.17,0.28,0.15,-0.03,-0.15,-0.07,0.059,0.08,0.032,-0.015,-0.02]';
f=fittype('a*cos(k*x)*exp(w*x)','independent','x');
cfun=fit(x,y,f)
xi=0:0.1:20;
yi=cfun(xi);
plot(x,y,'r*',xi,yi,'b-')

第四节插值问题

一、一元函数插值

（1）已知两个节点时，得线性插值多项式:

（2）已知三个节点时，得抛物插值多项式:

已知n+1个节点时，可得n次拉格朗日插值多项式。

若F(x)为分段函数，称为分段插值；插值中使用较多的是分段线性插值和三次样条插值。

分段线性插值

设函数y=f(x)的插值节点为(xi,yi)， (i=0,1,2, …,n) 将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数记为：In(x), 它满足In(xi)=yi ,且在每个小区间上是一次函数。

分段线性插值具有收敛性：但在节点处In(x)不可导，光滑性较差。

Matlab中的一维插值函数interp1

格式1：yi=interp1(x,y,xi) 功能：已知函数y=f(x)的一组插值节点(x,y)，求出其插值函数F(x)在点xi处的函数值yi。

格式2：yi=interp1(y,xi) 默认x=1:n,n为向量y的元素个数。

格式3：yi=interp1(x,y,xi,’method’)

Method：

(1)nearest 线性最近项插值

(2)linear 线性插值

(3)spline 三次样条插值

(4)cubic 立方插值

缺省时默认为为分段线性插值

例1 通过实验测得某函数的一组数据如下： x: 0,3,5,7,9,11,12,13,14,15 y: 0,1.2,1.65,2.1,2.15,2.0, 1.85,1.65,1.55,1.25 试作出其插值函数的图形

x=[0,3,5,7,9,11,12,13,14,15];

y=[0,1.2,1.65,2.1,2.15,2.0,1.85,1.65,1.55,1.25];

xi=0:0.1:15;

yi=interp1(x,y,xi,'spline');

yi1=interp1(x,y,xi,'linear');

yi2=interp1(x,y,xi,'cubic');

plot(x,y,'*',xi,yi,'r-',xi,yi1,‘b-',xi,yi2,‘g-')

legend('节点','三次样条插值','线性插值','立方插值')

三次样条插值是解决一维插值问题最常用的方法， Matlab中实现三次样条插值的方法有：

（1）yi=interp1(x,y,xi,’spline’)

（2）使用spline函数：

用法1：yi=spline(x,y,xi)（缺点：找不到方程表达式）效果同（1）

用法2：pp=spline(x,y)

获得三次样条插值的分段多项式pp，可使用ppval计算插值。

（3）使用csape函数：

用法： pp=csape(x,y) 可以添加参数选择边界条件。

例2 通过实验测得某函数的一组数据如下： x: 0,3,5,7,9,11,12,13,14,15 y: 0,1.2,1.65,2.1,2.15,2.0,1.85,1.65,1.55,1.25 试作出其插值函数的图形。

x=[0,3,5,7,9,11,12,13,14,15];

y=[0,1.2,1.65,2.1,2.15,2.0,1.85,1.65,1.55,1.25];

xi=0:0.1:15;

yi=interp1(x,y,xi,'spline');

yi1=interp1(x,y,xi,'linear');

yi2=interp1(x,y,xi,'cubic');

plot(x,y,'*',xi,yi,'r-',xi,yi1,‘b-',xi,yi2,‘g-') legend('节点','三次样条插值','线性插值','立方插值')

clear,clc

x=[0,3,5,7,9,11,12,13,14,15];

y=[0,1.2,1.65,2.1,2.15,2.0,1.85,1.65,1.55,1.25];

xi=0:0.1:15;

S=csape(x,y)

%Sa=spline(x,y)

P=S.coefs

% Pa=Sa.coefs（系数信息）

yi=ppval(S,xi);

plot(x,y,'o',xi,yi,'r');

二、二元函数插值

1. 网格节点数据插值函数：interp2

格式：z=interp2(x0,y0,z0,x,y,’method’)

x0,y0,z0（横纵竖）:插值节点坐标，要求x0,y0单调；

x,y是被插值点的横坐标与纵坐标( x,y不能超过x0,y0的范围)，z是被插值点的函数值。

Method：(1)nearest 最邻近插值 (2)linear 双线性插值 (3)cubic双三次插值默认为双线性插值。

4800

1350

1370

1390

1400

1410

960

940

880

800

690

570

430

290

210

150

4400

1370

1390

1410

1430

1440

1140

1110

1050

950

820

690

540

380

300

210

4000

1380

1410

1430

1450

1470

1320

1280

1200

1080

940

780

620

450

370

350

3600

1420

1430

1450

1480

1500

1550

1510

1430

1300

1200

980

850

750

550

500

3200

1430

1450

1460

1500

1550

1600

1550

1600

1550

1500

1550

1500

2800

950

1190

1370

1500

1200

1100

1550

1600

1550

1380

1070

900

1050

1150

1200

2400

910

1090

1270

1500

1200

1100

1350

1450

1200

1150

1010

880

1000

1050

1100

2000

880

1060

1230

1390

1500

1400

900

1100

1060

950

870

900

930

950

1600

830

980

1180

1320

1450

1420

1400

1300

700

900

850

840

380

780

750

1200

740

880

1080

1130

1250

1280

1230

1040

900

500

700

780

750

650

550

800

650

760

880

970

1020

1050

1200

830

800

700

300

500

550

480

350

400

510

620

730

800

850

870

850

780

720

650

500

200

300

350

320

370

470

550

600

670

690

670

620

580

450

400

300

100

150

250

y/x

400

800

1200

1600

2000

2400

2800

3200

3600

4000

4400

4800

5200

5600

clear,clc
 load data404
 x=0:400:5600;  y=4800:-400:0; 
 [X,Y]=meshgrid(x,y);
 surf(X,Y,Z)（原始数据的图形）
 pause
 figure(2)
 xi=linspace(0,5600,80);（取80个节点）
 yi=linspace(0,4800,80); 
 [Xi,Yi]=meshgrid(xi,yi);
 Zi=interp2(X,Y,Z,Xi,Yi); 
 surf(Xi,Yi,Zi)
 pause 
figure(3)
 t=0:100:1600; （画等高线）
 [c,h]=contourf(Xi,Yi,Zi,t);
 clabel(c,h)（标上等高线）
 colormap cool（加色调）
 colorbar（加色标）

2. 散点数据插值函数：griddata 已知n个插值节点(xi,yi,zi), (i=1,2,…,n), 求在点(x,y)处的插值z， matlab提供函数griddata。

格式：cz=griddata(x,y,z,cx,cy,’method’)

其中x,y,z 均为n 维向量，指明所给数据点（插值节点）的横坐标、纵坐标和竖坐标。 cx,cy是给定被插值点的横坐标和纵坐标，cz为相应点的竖坐标。

若cx,cy是向量，则给定以它们所确定网格点的横坐标和纵坐标，这时要求cx,cy一个为行向量一个为列向量。编程时也可先用meshgrid将cx,cy定义成网格矩阵。

clear ,clc
x=[129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5]; 
 y=[7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5];
 z=-[4,8,6,8,6,8,8,9,9,8,8,9,4,9];
 [x1,y1]=meshgrid(75:5:200,-50:5:150); 
 z1=griddata(x,y,z,x1,y1,’v4’);
 surf(x1,y1,z1)
 figure(2)
 [c,h]=contourf(x1,y1,z1);
 clabel(c,h)

第五节数值积分

一、数值积分的常用方法

（1）用梯形法计算积分，适用于被积函数为离散数据时，求函数的定积分。该函数调用格式：I=trapz(x,y)(分点越密效果越好)

求积分：

clc,clear
format long
ac=@(x)sin(x)./x
x1=pi/4:pi/50:pi/2;
y1=ac(x1);
s1=trapz(x1,y1)
x2=pi/4:pi/100:pi/2;
y2=ac(x2);
s2=trapz(x2,y2)

(2) 基于变步长辛普森法计算积分(了解) quad或quadtx，

调用格式： [I,n]=quad(‘fname’,a,b,Tol,trace)

其中：‘fname’是被积函数名。

a,b是积分上下限。

Tol是精度控制值，省却时取0.001

Trace:控制是否显示展现积分过程，取0不展现

I：积分值（返回量）

n：被积函数调用次数

如：求积分： ac=@(x)sin(x)./x s=quad(ac,pi/4,pi/2)

如：求积分： ac=@(x)sin(x)./x s=quad(ac, realmin（最小正实数）,pi/2)

做出函数的图形：

x=1:0.05:20;
n=length(x);
f=@(t)sin(t)./t;
for i=1:n
    y(i)=quad(f,0,x(i));
end
plot(x,y)

(3) 高精度Lobatto积分法，格式：z = quadl(Fun,a,b)

(4) 自适应Gauss-Kronrod数值积分法 z = quadgk(Fun,a,b)

如：求积分：

f=@(x)exp(-x.^2)
 y1=quadl(f,-1,1)
 y2=quadgk(f,-1,1)

(5) 积分法矢量化自适应simpson数值积分格式：z = quadv(Fun,a,b) 一次可以计算多个积分

如：求积分： F=@(x,n)1./((1:n)+x.^2); quadv(@(x)F(x,6),0,1)

含参变量的积分：做出函数

clc,clear
f=@(x,t)sin((x-t).^2);
t=-5:0.05:5;
n=length(t);
for i=1:n
s(i)=quad(f,0,pi,[],[],t(i));
end
plot(t,s)

mmk27

关注

0
点赞
踩
22

收藏

觉得还不错? 一键收藏
0
评论
matlab数值计算

第一节多项式运算一、多项式的表示、求值、求根 1. 表示： matlab中把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。 p=[2,3,-1,4,5]表示多项式：显示多项式的数学形式： p1=poly2str(p,'x') 如：p=[1,2,3,4] p1=poly2str(p,'x') 结果：p1 = x^3 + 2 x^2...
复制链接

扫一扫

专栏目录