昆明理工大学上机安排2数值计算.docx
《昆明理工大学上机安排2数值计算.docx》由会员分享，可在线阅读，更多相关《昆明理工大学上机安排2数值计算.docx（23页珍藏版）》请在冰点文库上搜索。
昆明理工大学上机安排2数值计算
内容：
本次上机主要练习常见的数值计算问题，包括线性代数、函数分析、数值微积分、常微分方程等，重点练习如何利用MATLAB提供的函数来实现数值计算，对于数学理论问题不做详细阐述，不清楚的请看相关数学书。
目的：
能应用MATLAB提供的强大函数进行复杂的方程组、微分、积分等运算。
注意：
MATLAB数值计算的结果为数值型数据，而不是数学上的解析表达式。
线性代数：
1.LU分解
一个矩阵可以分解为一个上三角矩阵和一个下三角矩阵的乘积，称之为LU分解。
LU分解是用高斯主元消去法实现的，通常要对主元位置进行交换，主元交换的方法是将被分解矩阵左乘一个由0-1构成的行交换阵。
【调用格式】
[L,U,P]=lu（X）对矩阵X进行LU分解，并进行主元交换，
[L,U]=lu（X）对矩阵X进行LU分解，无主元交换
【说明】L为主对角元素为1的下三角矩阵，U为上三角矩阵，P为行交换矩阵。
2.行列式和求逆
矩阵的行列式和求逆可以通过LU分解的方法求解。
【调用格式】
d=det（X）求矩阵X的行列式
Y=inv（X）求矩阵X的逆矩阵
例1对矩阵进行LU分解、求行列式和求逆操作，请查看计算结果。
>>A=[1,2,3;2,2,3;9,7,5];
>>[L1,U1]=lu（A）；%不带主元交换的LU分解，L1通常不是下三角阵
>>[L2,U2,P]=lu（A）%带主元交换的LU分解，L2为下三角阵
>>det（A）%行列式值
>>Y=inv（A）%矩阵求逆
3.特征值和特征向量
eig函数用于求解矩阵的特征值和特征向量。
【调用格式】
D=eig（A）计算矩阵A的特征值，D为特征值构成的向量
[V,D]=eig（A）计算矩阵A的特征值对角阵D和特征向量矩阵V
[V,D]=eig（A,'nobalance'）当矩阵A中有与截断误差近似的数值，用本指令
例2计算B的特征值及特征向量，观察结果。
>>B=[-2,1,1;0,2,0;-4,1,3];
>>D=eig（B）;
>>[V1,D1]=eig（B,'nobalance'）;%计算误差为1.0e-014数量级
4.线性方程组的解
在线性方程组中，独立方程的个数等于独立未知参数的个数，称为恰定方程；独立方程的个数大于独立未知参数的个数，称为超定方程；独立方程的个数小于独立未知参数的个数，称为欠定方程。
形如Ax=b的线性方程组可以用以下方法求解：
a.左除运算符法
左除运算符法的形式为:
x=A\b
对于一般的非奇异矩阵A，可以求得唯一数值解。
欠定方程和超定方程，可以获得最小二乘解。
b.广义逆法
如果用左除运算符求解的时候出现提示矩阵A为非奇异的警告或者解中出现Nan，则可以采用广义逆法，形式为:
x=pinv（A）*b
c.符号计算法
可以求得方程组的符号解，对于欠定方程可以求得具有自由变量的解。
例3求以下3个方程组的解,观察解是否正确。
I：
II：
III：
解：
方程I为恰定方程，用左除运算符可以求解
>>A1=[4,6,3;2,3,4;5,2,3];b1=[13;9;10];
>>x1=A1\b1
方程II为欠定方程，用3种方法分别求解
>>A2=A1;b2=b1;
>>A2（3,:
）=[];b2（3）=[];
>>x2=A2\b2%左除运算符法，求得一组特解，解中非零元素最少
x2=
0
1.6667
1.0000
>>x2=pinv（A2）*b2%广义逆法，求得一组特解
x2=
0.7692
1.1538
1.0000
>>x2=solve（'4*x+6*y+3*z=13,2*x+3*y+4*z=9','x,y,z'）;%符号运算法求通解
>>x2=[x2.x;x2.y;x2.z]
x2=
5/2-3/2*y
y
1
方程III为超定方程，用左除运算符法求最小二乘解
>>A3=A1;b3=b1;
>>A3（:
3）=[];
>>x3=A3\b3
x3=
1.4545
1.3636
函数分析
1.多项式的根
通过roots函数来求取多项式全部的根。
【调用格式】
r=roots（p）多项式求根函数
【说明】p为多项式的系数行向量，r为多项式所有根构成的列向量。
2.一元函数零点
fzero函数求取一元函数的精确零点。
【调用格式】
[x,fval,exitflag,output]=fzero（fun,x0,options）一元函数零点的完整调用格式
x=fzero（fun,x0）一元函数零点的最简调用格式
【说明】
1.fzero只能求得x0附近的单个零点，不能求取函数的所有零点。
2.输入变量fun表示一元函数，可以是字符串、内联函数或者函数句柄。
3.输入变量x0为零点的初始猜测值（自变量值）。
如果x0为标量，则求距离x0最近的那个零点；如果x0=[a,b]，要求fun（a）和fun（b）异号，此时求自变量在[a,b]区间内的零点。
4.输入变量options为优化迭代选项，是一个结构体。
5.输出变量x为零点处的自变量值，输出变量fval为零点处的函数值。
6.输出变量exitflag表示函数中止计算的条件。
若exitflag>0表示找到零点后退出。
输出变量output表示程序运行信息，是一个结构体。
例1求函数
在t>0区间内的所有零点。
解：
运行下列程序，请观察执行结果。
>>y=inline（'1-12.5*exp（-t）*sin（2*t+3.4）','t'）;%构造内联函数
>>t=0:
0.1:
10;%定义自变量取值区间
>>yv=feval（vectorize（y）,t）;%计算自变量区间的函数值
>>plot（t,yv）;
>>axistight;
>>gridon;
>>[tt,yy]=ginput%请在在函数曲线上用鼠标拾取所有零点，按回车键结束
tt=
yy=
0.0085
0.0085
>>[t1,fv1,flag]=fzero（y,tt
（1））;%求第一个拾取点的精确零点
>>[t2,fv2,flag]=fzero（y,tt
（2））;%求第二个拾取点的精确零点
附加学习及练习：
内联函数
用户可以用M文件来建立函数，函数的功能可以很复杂，函数的输出变量也可以有多个。
对于简单的数学表达式，用M文件来建立函数就显得不够方便。
MATLAB提供了内联函数的功能，内联函数可以将表达式转换为函数。
内联函数是MATLAB面向对象的一个类，其类型名为inline。
1.内联函数的建立
（1）g=inline（'expr'）将串表达式expr转换为内联函数
（2）g=inline（'expr','arg1','arg2',...）将串表达式expr转化为以arg1、arg2等为自变量（输入变量）的内联函数
（3）g=inline（'expr',n）将串表达式expr转化为以自变量x，P1，P2，…，Pn为自变量的内联函数。
其中P必须大写。
例阅读下列程序代码，了解inline的使用方法。
>>g=inline（'sin（alpha*x）','x','alpha'）
g=
Inlinefunction:
g（x,alpha）=sin（alpha*x）
>>g1=inline（'a*sin（x
（1））*cos（x
（2））','a','x'）%自变量为向量，函数值为标量
g1=
Inlinefunction:
g1（a,x）=a*sin（x
（1））*cos（x
（2））
>>g1（1,[pi/4,pi/4]）
ans=
0.5000
%自变量为向量，函数值为向量
>>g3=inline（'[sin（x
（1））;cos（x
（2））;sqrt（x
（1）.^2+x
（2）.^2）]'）
g3=
Inlinefunction:
g3（x）=[sin（x
（1））;cos（x
（2））;sqrt（x
（1）.^2+x
（2）.^2）]
>>g3（[3,4]）
ans=
0.1411
-0.6536
5.0000
>>g4=inline（'sqrt（x^2+P1^2+P2^2）',2）%inline的第三种格式的使用
g4=
Inlinefunction:
g4（x,P1,P2）=sqrt（x^2+P1^2+P2^2）
>>g4（1,3,5）
ans=
5.9161
>>feval（g4,1,3,5）%内联函数可以通过feval计算求值
ans=
5.9161
2.和内联函数有关的函数
class（fun）%获取内联函数的数据类型
char（fun）%获取内联函数的计算公式字符串
argnames（fun）%获取内联函数的输入变量名字
vectorize（fun）%使内联函数具有数组运算规则
字符串的求值
MATLAB提供了字符串求值的函数，利用这些函数，可以用字符串构造MATLAB的函数和命令，并运行这些字符串命令。
1.字符串表达式计算
【调用格式】
y=eval（'expression'）计算字符串表达式expression
[a1,a2,...]=eval（'function（b1,b2,...）'）计算函数调用的字符串表达式
【说明】eval的输入变量只能是字符串
例表达式字符串的计算。
解：
运行下列语句，观察执行结果。
clear
x=0.4;
cmd=['y=x*eye（3）,z=sin（x/2）+cos（x）'];
eval（cmd）;
who
2.字符串函数计算
【调用格式】
[y1,y2,...]=feval（'function',x1,...,xn）
【说明】
（1）'function'只能是函数名，不能是表达式字符串。
（2）x1、x2等是调用函数'function'的输入变量，即函数的自变量值。
（3）y1、y2等是函数的输出变量，即函数的返回值。
例feval的使用方法。
解：
执行下列语句，观察程序的运行结果
x=0:
2*pi/50:
2*pi;
feval（'plot',x,sin（x）,'r',x,cos（x）,'b:
'）;
%相当于运行plot（x,sin（x）,'r',x,cos（x）,'b:
'）
附加学习结束
3.多元函数的零点
fsolve函数求取多元函数的精确零点，但是必须提供零点的大致位置才能进行数值搜索。
【调用格式】
[x,fval,exitflag,output]=fsolve（fun,x0,options）多元非线性方程求解
x=fsolve（fun,x0）多元非线性方程求解的最简格式
【说明】
1.输入变量fun表示自变量为向量的函数，可以是字符串、内联函数或者函数句柄。
2.输入变量x0为零点的初始猜测向量（自变量值）。
3.输出变量x为零点处的自变量值，输出变量fval为零点处的函数值。
例求方程组
在x1=-5,x2=-5附近的解。
解：
（1）采用字符串表达函数
>>fun='[2*x
（1）-x
（2）-exp（-x
（1））,-x
（1）+2*x
（2）-exp（-x
（2））]';%字符串方式表达函数
>>x0=[-5;5];
>>[x,fval]=fsolve（fun,x0）
（2）采用内联函数表达函数
>>fun=inline（'[2*x
（1）-x
（2）-exp（-x
（1））,-x
（1）+2*x
（2）-exp（-x
（2））]','x'）;
>>x0=[-5;5];
>>[x,fval]=fsolve（fun,x0）
4.函数的极值点
MATLAB提供了3个求极值点的函数，其输入输出参数的定义和fsolve函数基本相同。
【调用格式】
[x,fval,exitflag,output]=fminbnd（fun,x1,x2,options）
【说明】求一元函数fun在自变量（x1，x2）区间的最小值
【调用格式】
[x,fval,exitflag,output]=fminsearch（fun,x0,options）
【说明】：
用单纯形法求多元函数fun在自变量向量x0附近的极小值点
【调用格式】
[x,fval,exitflag,output]=fminunc（fun,x0,options）
【说明】：
用拟牛顿法求多元函数fun在自变量向量x0附近的极小值点
例求二元函数
在x=-1.2,y=1附近的极小值。
解：
运行下列程序，观察运行结果。
>>fun=inline（'100*（x
（2）-x
（1）^2）^2+（1-x
（1））^2','x'）;
>>[x,fval]=fminsearch（fun,[-1.2,1]）
x=
1.00001.0000
fval=
8.1777e-010
5.数值微分
数值导数和微分是基于差分定义的，对原始数据的采样间隔依赖很大，如果原始数据含有噪声，则数值导数结果没有什么价值，因此尽量避免求数值导数，diff函数用于求相邻数据间的数值差分。
【调用格式】
DX=diff（X）求X相邻元素的一阶差分，即用后一个元素值减去当前元素值
DX=diff（X,n）求X相邻元素的n阶差分
DX=diff（X,n,dim）在dim指定的维上，求X相邻元素的n阶差分
【说明】
1.如果X是向量，差分是相邻元素相减；如果X是矩阵，差分是相邻行间对应元素相减。
2.差分数据DX元素的个数要比被操作数X少一个。
3.DX只是差分数据，如果相邻数据点之间的步长为DH，则DX./DH为导数数据。
例绘制函数f（x）=sinx的导函数在（0,2π）区间的曲线。
解：
运行下列程序，观察所得图形。
>>dx=2*pi/50;x=0:
dx:
2*pi;y=sin（x）;
>>dy=diff（y）;%求差分数据
>>df=dy./dx;%求导数数据
>>y=y（1:
end-1）;x=x（1:
end-1）;%差分数据少1个，裁剪原始函数的相关数据
>>plot（x,y,'r--',x,df,'b'）,axistight,gridon
6.数值积分
数值积分分为开型积分和闭型积分，二者的区别在于是否计算积分区间端点处的函数值。
MATLAB提供的数值积分函数有些适用于开型积分和闭型积分，有些只能用于闭型积分。
积分计算可以采用本节介绍的相关函数，也可以采用样条积分法，还可以采用符号积分法，后两种方法在下次上机介绍。
a.一元函数的数值积分
求一元函数数值积分的函数有很多，每个函数采用不同的数值算法，各有优缺点，精度和运算速度也不尽相同，quadl函数为最常用的数值积分函数。
【调用格式】
q=quadl（fun,a,b,tol,trace）采用递推自适应Lobatto法计算一元函数的积分
【说明】
（1）fun为被积函数，可以用字符串、内联函数和函数M文件的函数句柄表示，且被积函数表达式要按照数组运算规则来编写。
通常积分变量采用字母x。
（2）a和b为积分变量的积分上下限，为常数数值。
（3）tol为绝对误差，是一个标量，可以省略。
（4）trace为跟踪标志，当trace为非零值时，随积分进程会逐点画出被积函数。
（5）返回值q为数值积分的结果。
例求积分
解：
运行下列程序，观察运行结果。
>>fun=inline（'1./（x.^3-2.*x+5）','x'）;
>>q=quadl（fun,0,2）
b.二重数值积分
形如
的闭型二重数值积分的函数为dblquad。
【调用格式】
q=dblquad（fun,x1,x2,y1,y2,tol,method）
【说明】
（1）fun为被积函数f（x,y），可以用字符串、内联函数和函数M文件的函数句柄表示。
积分变量用x，y表示，x为内重积分变量，y为外重积分变量。
（2）x2和x1是变量x的积分上下限，y2和y1是变量y的积分上下限。
（3）method表示选用的积分算法，可以省略，缺省算法是@quad，还可以选则@quadl或者用户自定义的积分算法M函数文件的函数句柄。
（4）本函数无法计算内重积分上下限为函数表达式的情况。
c.三重数值积分
求形如
的闭型三重数值积分的函数为triplequad。
【调用格式】
q=triplequad（fun,x1,x2,y1,y2,z1,z2,tol,method））
【说明】
（1）fun为被积函数f（x,y,z），可以用字符串、内联函数和函数M文件的函数句柄表示。
积分变量用x，y，z表示，从内到外积分变量依次为x，y，z。
（2）x2和x1是变量x的积分上下限，y2和y1是变量y的积分上下限，z2和z1是变量z的积分上下限。
（3）其他事项同函数dblquad。
例计算定积分
解：
运行下列程序，观察运行结果。
>>fun=inline（'y*sin（x）+z*cos（x）','x','y','z'）;
>>q=triplequad（fun,0,pi,0,1,-1,1）
数据拟合
已知某些离散的原始数据，建立一个曲线方程，让它以最佳的方式反映原始数据的变化趋势，尽量避免出现局部的波动，这种方法称为拟合。
不要求拟合曲线经过所有的原始数据，最佳方式通常是指用拟合得到的数学模型计算出来的计算数据和原始数据之间的误差的平方和最小，这种方式称为最小二乘。
通常是已知曲线方程的形式（数学模型的结构），根据原始数据来计算曲线方程的参数（数学模型的参数）。
1.多项式拟合
多项式拟合是用一个n阶多项式模型去拟合原始数据的方法，需要计算的模型参数包括多项式的阶次和多项式的系数，polyfit函数实现多项式拟合。
【调用格式】
p=polyfit（x,y,n）根据给定的数据（x,y），计算n阶拟合多项式的系数向量p
ye=polyval（p,x）计算自变量为x时多项式p的值（估计值）
【说明】多项式拟合一般不要超过5阶，否则计算误差变大；拟合多项式只在原始数据范围内可以保证精度，超出范围使用拟合多项式无法保证预报的精度。
例已知五个数据点：
[1,5.5],[2,43],[3,128],[4,290],[5,498]，试画出这五个点拟合的三次曲线。
解：
运行下列程序，请观察图形特点。
>>x=[1,2,3,4,5];y=[5.5,43,128,290,498];
>>p=polyfit（x,y,3）
>>x2=1:
0.1:
5;y2=polyval（p,x2）;
>>plot（x,y,'o',x2,y2）,gridon
2.最小二乘拟合
a.线性最小二乘拟合
对于线性数学模型的参数估计，可以用形如Y=Ax+b的一阶多项式拟合来估计参数。
某些非线性模型经过变量替换也可以转换为线性模型，也可以采用线性估计方法，lsqlin函数实现线性最小二乘估计。
【调用格式】
a=lsqlin（Fun,a0,lb,ub,options,p1,p2,...）
b.非线性最小二乘估计
可以使用MATLAB提供的lsqnonlin函数实现非线性最小二乘估计。
【调用格式】
[a,resnorm,residual,exitflag,output]=lsqnonlin（Fun,a0,lb,ub,options,p1,p2,...）
a=lsqnonlin（Fun,a0）
【说明】
（1）函数的功能是求取向量的值，使向量函数
满足误差平方和最小
（2）输入参数Fun是被解函数，Fun是向量函数，其自变量a是向量。
Fun可以用字符串、内联函数或者函数M文件的函数句柄表示。
（3）输入变量a0为向量，表示所求变量a的初始猜测值。
（4）输入变量lb表示所求变量a的下限，输入变量ub表示所求变量a的上限。
如果没有上下限限制，可以用空矩阵[]表示，如果没有下限，可以用-inf表示。
（5）从第六个输入变量开始的参数p1,p2等是向Fun函数传递的参数，其名字应该和定义Fun函数的输入变量名相一致。
（6）输出变量a表示所求变量的值。
（7）如果将构造为残差函数，且为被估计参数向量，则lsqnonlin就变为求解非线性最小二乘估计的函数。
例混凝土的抗压强度随养护时间的延长而增加，其数学模型为
现将一批混凝土作成12个测试块，记录了养护时间x（日）及抗压强度y（kg/cm2）的数据，如表1所示。
试估计该数学模型的参数a1，a2，a3，a4。
表1混凝土的抗压强度y和养护时间x的实测数据
x
2
3
4
5
7
9
12
14
17
21
28
56
y
35
42
47
53
59
65
68
73
76
82
86
99
解：
（1）编写残差计算函数文件Residuals.m
functionE=Residuals（a,x,y）
x=x（:
）;y=y（:
）;%将原始数据变为列向量
Y=a
（1）+a
（2）*exp（a（4）*x）+a（3）*exp（-a（4）*x）;%估计参数下的计算值
E=Y-y;%计算值减去原始值等于残差
（2）估计模型参数，并绘制模型曲线，如图5.3.2所示
>>x=[234579121417212856];
>>y=[354247535965687376828699];
>>a0=[50,0,-50,0];
>>options=optimset（'lsqnonlin'）;
>>options.TolX=0.01;options.Display='off';
>>a=lsqnonlin（@Residuals,a0,[],[],options,x,y）
a=
92.36270.0500-65.28500.0878
>>xx=min（x）:
max（x）;
>>yy=a
（1）+a
（2）*exp（a（4）*xx）+a（3）*exp（-a（4）*xx）;%估计参数下的计算值
>>plot（x,y,'o',xx,yy,'r'）
3.插值和样条
曲线拟合是用带
升级会员 