版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
科学计算与MATLAB主讲:唐建国中南大学材料科学与工程学院2013.11第一章绪论1科学计算数值计算的定义数值计算的研究过程数值计算的误差分析数值计算的稳定性与收敛性2MATLAB简介MATLAB的类别MATLAB的概述3课程的要求与学习方法课程的内容课程的要求课程的学习与考核教材与主要参考书小结如果存在一适当小的正数εr
,使得则称εr为相对误差限。
例:x=15, ε(x)=2;εr(x)=2/15=13.33%
y=1000,ε(y)=5,εr(y)=5/1000=0.5%2.有效数字:定义:若近似值x的误差限是某一位的半个单位,该位到x的第一位非零数字共有n位,就说x有n位有效数字,x可表示为其中,a1≠0.且a1,a2,…,an
都是0~9中的任一整数。其绝对误差限满足:例.
测量一物体的长度为954cm,问测量数据的相对误差限多大?解:因实际问题所截取的近似数,其绝对误差限一般不超过最小刻度的半个单位,故当x=954cm时,有ε(x)=0.5cm, 而x的相对误差
er(x)≤0.5/954=0.0005241….<0.00053=0.053%
故εr(x)=0.053%.第二讲MATLAB基础知识(数值、符号计算)1MATLAB的启动与运行Matlab的启动和退出Matlab的工作窗口Matlab的指令行的操作Matlab的帮助系统2MATLAB的矩阵与数值计算Matlab数据类型Matlab的变量及其命名Matlab矩阵及其运算Matlab的数组关系/逻辑运算Matlab的多项式运算3MATLAB的符号计算小结1.5Matlab的帮助系统MATLAB帮助命令包括help、lookfor以及模糊查询1help命令在命令窗口中输入help命令将会显示当前帮助系统中所包含的所有项目,即搜索路径中所有的目录名称。也可以通过help加函数名来显示该函数的帮助说明。2lookfor命令对搜索范围内的M文件的第一行进行关键字搜索,条件比较宽松。若对M文件进行全文搜索,加上-all选项。3模糊查询
MATLAB6.0以上的版本提供了一种类似模糊查询的命令查询方法,只需要输入命令的前几个字母,然后按Tab键,系统就会列出所有以这几个字母开头的命令。
数值显示格式任何MATLAB的语句的执行结果都可以在屏幕上显示,同时赋值给指定的变量,没有指定变量时,赋值给一个特殊的变量ans,数据的显示格式由format命令控制。
format只是影响结果的显示,不影响其计算与存储;MATLAB总是以双字长浮点数(双精度)来执行所有的运算。
如果结果为整数,则显示没有小数;如果结果不是整数,则输出形式有:format(short):短格式(5位定点数)99.1253formatlong:长格式(15位定点数)
99.12345678900000formatshorte:短格式e方式9.9123e+001formatlonge:长格式e方式
9.912345678900000e+001formatbank:2位十进制99.12formathex:十六进制格式符号变量和符号表达式的创建
MATLAB提供了两个建立符号量的命令:syms
和sym
说明:a1a2a3...为需要定义为符号量的标识符,不能是数字、函数表达式或方程式。a1a2a3...均不用引号界定,而用空格分开。1用syms创建符号量格式:symsa1a2a3…说明:输入参量A可以是数字、字符串也可以是字符串变量名、字符表达式或字符方程。回车后将A定义成符号量并赋值给标识符。2用sym创建符号量、符号表达式
格式:标识符=sym(A)3.矩阵算法和数组算法矩阵算法
把矩阵看作一个整体,各种运算完全按照线性代数代表的矩阵运算法则进行,运算的书写形式和运算符号都与矩阵理论完全相同。数组算法把矩阵看作由其元素构成的一组数据(数组),各种运算是在参与运算矩阵的对应元素之间进行的数与数的运算,这种运算方便对大批数据的处理和一次求出多个函数值。数组算法的运算符主要有.*./.\.^
Matlab运算符
算术运算符操作符说明+加A+BAB必须大小相同,或一个是标量-减A-BAB必须大小相同,或一个是标量*矩阵乘A*BA的列数等于B的行数,一个可以是标量.*数组乘A.*BAB必须大小相同,一个可以是标量\矩阵左除A\B=A-1*B等效于A*X=B求Xinv(A).\数组左除A.\BBij/Aij/矩阵右除A/B=A*B-1等效于X*B=A求X./数组右除A./BAij/Bij^矩阵乘方A^pA自乘p次.^数组乘方A.^pA中每个元素的p次方如:a=[12;34];b=[35;59]》c=a+bd=a-b》c=d=47-2-3813-2-5》a*b=[1323;2951]》a/b=[-0.500.50;3.50–1.50]》a\b=[-1-1;23]》a^3=[3754;81118]》a.*b=[310;1536]》a./b=[0.330.40;0.600.44]》a.\b=[3.002.50;1.672.25]》a.^3=[18;2764]
只有维数相同的矩阵才能进行加减运算。注意只有当两个矩阵中前一个矩阵的列数和后一个矩阵的行数相同时,才可以进行乘法运算。a\b运算等效于求a*x=b的解;而a/b等效于求x*b=a的解。只有方阵才可以求幂。转置:对于实矩阵用(’)符号或(.’)求转置结果是一样的;然而对于含复数的矩阵,则(’)将同时对复数进行共轭处理,而(.’)则只是将其排列形式进行转置。4.常用矩阵运算逆矩阵与行列式计算求逆:inv(A);求行列式:det(A)要求矩阵必须为方阵5、矩阵分解(1)奇异值分解[U,S,V]=svd(A)例:a=9868可以验证:u*u’=Iv*v’=Iu*s*v’=a求矩阵A的奇异值及分解矩阵,满足U*S*V’=A,其中U、V矩阵为正交矩阵(U*U’=I),S矩阵为对角矩阵,它的对角元素即A矩阵的奇异值。[u,s,v]=svd(a)u=0.7705-0.63750.63750.7705s=15.5765001.5408v=0.6907-0.72310.72310.6907(2)特征值分解[V,D]=eig(A)例:
a=9868[v,d]=eig(a)v=0.7787-0.73200.62740.6813d=15.4462001.5538求矩阵A的特征向量V及特征值D,满足A*V=V*D。其中D的对角线元素为特征值,V的列为对应的特征向量。如果D=eig(A)则只返回特征值。(3)正交分解[Q,R]=qr(A)例:a=9868[q,r]=qr(a)q=-0.8321-0.5547-0.55470.8321r=-10.8167-11.094002.2188将矩阵A做正交化分解,使得Q*R=A,其中Q为正交矩阵(其范数为1,指令norm(Q)=1),R为对角化的上三角矩阵。(4)三角分解[L,U]=lu(A)将A做对角线分解,使得A=L*U,其中L为下三角矩阵,U为上三角矩阵。注意:L实际上是一个“心理上”的下三角矩阵,它事实上是一个置换矩阵P的逆矩阵与一个真正下三角矩阵L1(其对角线元素为1)的乘积。[L1,U1,P]=lu(A)例:a=[123;456;789]
比较:[l1,u1,p]=lu(a)[l,u]=lu(a)3、MATLAB的符号计算
A=sym(A)sum(A)sum(A’)’
det(A)inv(A)
eig(A)
svd(A)符号计算常规:第三讲MATLAB基础知识(图形、程序设计)1MATLAB的图形功能二维图形的基本函数三维图形的基本函数图形窗口函数简单动画函数2MATLAB的程序设计M文件介绍Matlab的程序结构Matlab的文件操作Matlab的函数操作小结1.1.1基本二维图形函数1:plot:绘制二维数据图形格式:plot(X,′S′)plot(X,Y,′S′)plot(X1,Y1,′S1′,X2,Y2,′S2′,......,X3,Y3,′S3′)说明:参数′S′控制数据点的标记、曲线类型和曲线色彩,三者置于一对单引号内。常用的绘图选项选项含义选项含义-实线*用星号标出数据点--虚线.用点号标出数据点
:点线o用圆圈号标出数据点-.点划线x用叉号标出数据点b蓝色+用加号标出数据点g绿色s用正方形标出数据点r红色D用菱形标出数据点c青色V用下三角标出数据点m洋红^用上三角标出数据点y黄色<用左三角标出数据点k黑色>用右三角标出数据点w白色H用六角形标出数据点P用五角形标出数据点函数2:
fplot:绘制y=f(x)图形格式:fplot(fname,lims,′S′)说明:其中fname为函数名或单引号界定的函数表达式,lims为x,y的取值范围,′S′定义与plot函数相同。函数3:
ezplot:绘制隐函数图形格式:ezplot(f,[xmin,xmax,ymin,ymax])说明:在区间xmin<x<xmax和ymin<y<ymax绘制f(x,y)=0的图形。,默认区间-2π<x<2π和-2π<y<2π1f(x,y)=0格式:ezplot(X,Y,[tmin,tmax])说明:在区间tmin<t<tmax绘制x=X(t)和y=Y(t)的图形,默认区间0<t<2π2x=X(t)y=Y(t)格式:ezplot(f,[a,b])说明:在区间a<x<b绘制y=f(x)的图形,默认区间-2π<x<2π3y=f(x)2特殊坐标二维图形(1)极坐标曲线格式:polar(theta,rho,′S′)theta:角度向量rho:幅值向量′S′:控制参数theta=0:0.1:8*pi;r=cos(4*theta)+1/4;polar(theta,r)例(2)对数坐标曲线x=0:0.01:5;y=10.^x;plot(x,y),gridon函数功能semilogxx轴对数坐标,y轴线性坐标semilogyx轴线性坐标,y轴对数坐标loglogxy轴均为对数坐标例x=0:0.01:5;y=10.^x;semilogy(x,y)gridon3应用型绘图指令可用于数值统计分析或离散数据处理
bax(x,y);hist(y,x);stairs(x,y);stem(x,y)
例:x=1:10;y=rand(10,1);subplot(2,2,1)bar(x,y);subplot(2,2,2)hist(y,x);subplot(2,2,3)stairs(x,y);subplot(2,2,4)stem(x,y);1图形标注函数中的说明文字,除使用标准的ASCII字符外,还可使用LaTeX格式的控制字符,这样就可以在图形上添加希腊字母、数学符号及公式等内容。例如,text(0.3,0.5,′sin({\omega}t+{\beta})′)将得到标注效果sin(ωt+β)。上述函数除legend外,均可以用于三维函数。主要函数:title(′图形名称′)xlabel(′x轴说明′)ylabel(′y轴说明′)text(x,y,′图形说明′)legend(′图例1′,′图例2′,...)1.1.2二维图形的控制2坐标控制axis函数
主要格式axis([xminxmaxyminymaxzminzmax])axisequal:纵、横坐标轴采用等长刻度。axissquare:产生正方形坐标系(缺省为矩形)。axisauto:使用缺省设置。axisoff:取消坐标轴。axison:显示坐标轴。gridon/off:控制是否画网格线。boxon/off:控制是否加边框线。holdon/off
控制是否刷新当前轴及图形1三维曲线图t=-pi:0.1:8*pi;x=sin(t);y=cos(t);plot3(x,y,t,′-r′)xlabel(′sin(t)′);ylabel(′cos(t)′);zlabel(′t′)plot3函数格式:Plot3(x1,y1,z1,′S1′,x2,y2,z2,′S2′,…xn,yn,zn,′Sn′)
1.2三维图形的基本函数2三维网格图meshgrid函数:产生平面区域内的网格坐标矩阵。格式:[X,Y]=meshgrid(A,B)说明:语句执行后,矩阵X的每一行都是向量A,行数等于向量B的元素的个数,矩阵Y的每一列都是向量B,列数等于向量A的元素的个数。mesh函数格式:mesh(x,y,z)说明:一般情况下,x,y,z是维数相同的矩阵。[x,y]=meshgrid(0:0.08:2*pi);z=sin(x).*cos(y);figure(1)mesh(x,y,z)xlabel('x'),ylabel('y')zlabel('sin(x)cos(x)')gridon,boxonfigure(2)mesh(z),boxonsurf函数格式:surf(x,y,z)[x,y]=meshgrid(0:0.08:2*pi);z=sin(x).*cos(y);surf(x,y,z)zlabel(′sin(x)cos(x)′)gridonboxon1.3图形窗口函数在实际应用中,有时需要在不同图形窗口或一个图形窗口中绘制若干个独立的图形,这就需要选取不同的图形窗口或对图形窗口分割。函数1:figure格式:figure(n)说明:该函数打开不同的图形窗口。n为图形窗口排序号。默认时打开的是1号图形窗,即当前窗。函数2:subplot格式:subplot(m,n,p)说明:该函数将当前图形窗口分成m×n个绘图区,即每行n个,共m行。区号按行优先编号,且选定第p个区为当前活动区。在每一个绘图区允许以不同的坐标系单独绘制图形。2.2Matlab的程序结构matlab的程序结构与c语言相似,主要有以下三种基本结构顺序结构选择结构循环结构2选择结构(1)if语句a单分支if语句:
格式
if条件语句组
endb双分支if语句:
格式
if条件语句组1else
语句组2endc多分支if语句
格式
if条件1
语句组1
elseif
条件2
语句组2……
elseif
条件m
语句组melse
语句组nend(2)switch语句switch表达式
case表达式1
语句组1case表达式2
语句组2……case表达式m
语句组motherwise
语句组nend说明:当表达式的值等于表达式1的值时,执行语句组1。当表达式的值等于表达式m的值时,执行语句组m。当表达式的值不等于case所列的表达式的值时,执行语句组n。当任意一个分支的语句执行完后,直接执行switch语句的下一句。(3)try语句try
语句组1catch
语句组2end说明:
try语句先试探性执行语句组1,如果语句组1在执行过程中出现错误,则将错误信息赋给保留的lasterr变量,并转去执行语句组2。3循环结构(1)for语句格式
for循环变量=表达式1:表达式2:表达式3
循环体语句
end说明:表达式1的值为循环变量的初值,表达式2的值为步长,表达式3的值为循环变量的终值。步长为1时,表达式2可以省略。格式
for循环变量=矩阵表达式循环体语句
end说明:执行过程是依次将矩阵的各列元素赋给循环变量,然后执行循环体语句,直至各列元素处理完毕,循环次数为列数。s=0;a=[123;456];fork=as=s+k;enddisp(s)(2)while语句格式
while(条件)
循环体语句
end说明:若条件成立,则执行循环体语句,执行后再判断条件是否成立,如果不成立则跳出循环。(3)break和continue语句break语句:用于终止循环的执行。当在循环体内执行到该语句时,程序将跳出循环,继续执行循环语句的下一语句。continue语句:控制跳过循环体中的某些语句。当在循环体内执行到该语句时,程序将跳过循环体中所有剩下的语句,继续下一次循环。2.4Matlab的函数操作1函数的调用说明:函数调用时各实参出现的顺序、个数,应与函数定义时形参的顺序、个数一致,否则会出错。函数调用时,先将实参传递给相应的形参,从而实现参数传递,然后再执行函数的功能。格式
[输出实参表]=函数名(输入实参表)2函数参数的可调性说明:在调用函数时,用nargin和nargout分别记录调用该函数时的输入实参和输出实参的个数。只要在自定义函数文件中包含这两个函数,就可以准确地知道该函数文件被调用时的输入输出参数个数,从而决定函数如何进行处理。
nargin
和nargout1、在动手编程之前,明确程序的目的,设想解决方案,作出初步的流程图。2、%后面的内容是程序的注解,要善于运用注解使程序更具可读性。3、在主程序开头用clear指令清除变量的习惯。但在子程序中不要用。4、参数值要集中放在程序的开始部分,以便维护。5、程序尽量模块化。6、充分利用Debugger来进行程序的调试(设置断点、单步执行、连续执行),并利用其他工具箱或图形用户界面(GUI)的设计技巧,将设计结果集成到一起。7、设置好MATLAB的工作路径,以便程序运行。8、input指令可以用来输入一些临时的数据;而对于大量参数,则通过建立一个存储参数的子程序,在主程序中用子程序的名称来调用。程序设计的注意事项
3.
请补充语句以画出如图所示的图形:[x,y]=meshgrid(-2:0.1:2,-2:0.1:2);Z=x.*exp(-x.^2-y.^2);
;a)Plot3(x,y,Z)b)plot3(x,y,Z)c)mesh(x,y,Z)d)plot3(x,y,z)第四讲插值法引言Lagrange插值分段低次插值Hermite插值三次样条插值插值方法比较MATLAB的插值函数小结7、MATLAB的插值函数插值函数及其功能函数功能描述Interp1一维插值Interp2二维插值Interp3三维插值Interpft一维快速傅立叶插值InterpnN-D维插值Spline
立体样条数据插值格式:yi=inter1(x,y,xi,′method′)说明:xy输入插值节点向量xiyi
插值点methodnearest:最近插值,用直角折线连接节点linear:线性插值,参数省略时,默认此项pehip
:分段三次Hermite插值,具有一阶导数连续splin:三次样条插值,具有一阶、二阶导数连续第五讲函数逼近与拟合法引言函数逼近傅里叶逼近最小二乘法拟合最小二乘法多元线性拟合非线性拟合MATLAB的拟合函数小结将其表示成矩阵形式:其系数矩阵为对称阵。所以正规方程组的系数矩阵非奇异,即根据Crame法则,正规方程组有唯一解,称其为最小二乘解。function[a,b]=LZXEC(x,y)%离散试验数据点的线性最小二乘拟合%试验数据点的x坐标向量:X%试验数据点的y坐标向量:Y%拟合的一次项系数:a%拟合的常数项:bif(length(x)==length(y))n=length(x);else
disp('x和y的维数不相等!');return;end%维数检查A=zeros(2,2);A(2,2)=n;B=zeros(2,1);fori=1:nA(1,1)=A(1,1)+x(i)*x(i);A(1,2)=A(1,2)+x(i);B(1,1)=B(1,1)+x(i)*y(i);B(2,1)=B(2,1)+y(i);endA(2,1)=A(1,2);s=A\B;a=s(1);b=s(2);例X12345y1.51.843.45.7>>x=1:5;>>y=[1.51.843.45.7];>>[a,b]=LZXEC(x,y)a=1.0000b=0.2800常用的线性变换函数变换后的函数4、matlab拟合函数线性拟合函数格式:p=linefit(x,y)说明:x,y输入同维数据向量p输出拟合多项式的系数向量多项式拟合函数格式:p=polyfit(x,y,m)说明:x,y输入同维数据向量m拟合多项式次数p输出拟合多项式的系数向量4.设y=span{1,x,x2},用最小二乘法拟合如下表数据。x0.51.01.52.02.53.0y1.752.453.814.808.008.60计算的结果为:a)0.49001.25010.8560b)0.85601.25010.4900c)-0.63413.8189-3.7749d)3.8189-3.77492.8533>>x=0.5:0.5:3.0;>>y=[1.75,2.45,3.81,4.80,8.00,8.60];>>a=polyfit(x,y,2)第六讲数值积分与微分引言数值积分矩形积分近似计算梯形积分近似计算抛物线形积分近似计算牛顿-科茨(Newton-Cotes)公式自适应(Simpson)求积法高斯(Gauss)求积法数值微分MATLAB的积分和微分函数小结MATLAB中提供的积分函数:矩形积分近似计算梯形积分近似计算抛物线形积分近似计算牛顿-科茨(Newton-Cotes)公式高斯(Gauss)求积法自适应(Simpson)求积法龙贝格积分法样条函数求积分法简单奇异积分重积分计算
4.2matlab数值计算定积分格式:Z=trapz(x,y)(1)梯形求积公式说明x,y分别代表数目相同的阵列或矩阵。y与x具有一定的函数关系。格式:[I,n]=quad('fname',a,b,tol,trace)(2)自适应辛卜生法说明返回参数I即定积分值,n为被积函数的调用次数。fname是被积函数名。定义被积函数须用数值运算符。a和b分别是定积分的下限和上限。tol用来控制积分精度,缺省时取tol=0.001。trace控制是否展现积分过程,若取非0则展现积分过程,取0则不展现,缺省时取trace=0。格式:[I,n]=quad8('fname',a,b,tol,trace)(1)Newton-Cotes法说明tol用来控制积分精度,缺省时取tol=10-6。其它与quad相同。该函数可以更精确地求出定积分的值,且计算精度相同条件下,函数调用的步数明显小于quad函数,从而保证能以更高的效率求出所需的定积分值。格式:[I,n]=quadl('fname',a,b,tol,trace)(4)高斯-洛巴托(Gauss-Labatto)法格式:I=dblquad('fname‘,xmin,xmax,ymin,ymax,tol,trace)(5)二重积分
4.3matlab数值计算微分格式:f=diff(fun)F=diff(fun,’x’)F=diff(fun,’x’,n)F=diff(S)F=diff(S,’x’)F=diff(S,’x’,n)(1)数值微分和差分一元函数求导格式:[gx,gy]=gradient(F)[gx,gy]=gradient(F,H)[nx,ny,nz]=surfnorm(X,Y,Z)二元、多元函数求导clearall;df1=diff('sin(x)');df2=diff('sqrt(x)')df3=diff('log(x)')df4=diff('log(x*y)','x')df5=diff('sin(x*y)','y')df6=diff('sin(x*y)','x',3)Ans:df1=cos(x)df2=1/2/x^(1/2)df3=1/xdf4=1/xdf5=cos(x*y)*xdf6=-cos(x*y)*y^3clearall;F=[11.21.42.35;0-0.6342;-177.291.4];[gx,gy]=gradient(F)%计算梯度n=surfnorm(F)%计算法向量functionf=fun(x)f=x.*x.*sqrt(2*x+3)quadl(‘fun’,1,3,1e-10)
或quadl('exp(-x/2)',1,3,1e-10)或fun=@(x)exp(-x/2);quadl(fun,1,3,1e-10)symsxf=(x^3+x^2+x+1)^(1/3)-sqrt(x^2+x+1)*log(exp(x)+x)/xlimit(f,x,inf)symsxy=log(1+x)f=diff(y,2)subs(f,1)第七讲常微分方程数值解法引言欧拉近似方法龙格-库塔(R-K)方法MATLAB的常微分方程函数小结4、MATLAB的常微分方程函数ode45ode23ode113ode15sode23sode23tode23tb格式[x,y]=ode45(′fun′,[x0,xn],y0,option]说明:适用于求解一阶常微分方程组fun定义微分方程组的函数文件名[x0,xn]求解区域y0初始条件向量option可选参数,由ODESET函数设置,比较复杂x输出自变量向量,y输出[y,y′,y″,..]没有一种算法可以有效地解决所有的ODE问题,因此MATLAB提供了多种ODE函数。函数ODE类型特点 说明ode45非刚性单步法;4,5阶R-K方法;累计截断误差为(△x)3大部分场合的首选方法ode23非刚性单步法;2,3阶R-K方法;累计截断误差为(△x)3使用于精度较低的情形ode113非刚性多步法;Adams算法;高低精度均可到10-3~10-6
计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法;Gear’s反向数值微分;精度中等若ode45失效时,可尝试使用ode23s刚性单步法;2阶Rosebrock
算法;低精度当精度较低时,计算时间比ode15s短ode23tb刚性梯形算法;低精度当精度较低时,计算时间比ode15s短例题:用MATLAB的符号解法,求解常微分方程:>>dsolve('Dy+3*x*y=x*exp(-x^2)')ans=(1/3*exp(-x*(x-3*t))+C1)*exp(-3*x*t)
>>dsolve('Dy+3*x*y=x*exp(-x^2)','x')ans=exp(-x^2)+exp(-3/2*x^2)*C1例题:采用ODE45求解描述某刚性问题的微分方程:functiondy=rigid(t,y)dy=zeros(3,1);%行向量dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);options=odeset('RelTol',1e-4,'AbsTol',[1e-41e-41e-5]);[T,Y]=ode45(@rigid,[012],[011],options);plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')legend('y1','y2','y3')例题:用MATLAB函数ode23,ode45,ode113,求解多阶常微分方程:functiondy=myfun03(x,y)dy=zeros(3,1)%初始化变量dydy(1)=y(2);%dy(1)表示y的一阶导数,其等于y的第二列值dy(2)=y(3);%dy(2)表示y的二阶导数dy(3)=2*y(3)/x^3+3*y(2)/x^3+3*exp(2*x)/x^3%dy(3)表示y的三阶导数%用ode23ode45ode113解多阶微分方程clear,clc[x23,y23]=ode23('myfun03',[1,10],[11030]);[x45,y45]=ode45('myfun03',[1,10],[11030]);[x113,y113]=ode113('myfun03',[1,10],[11030]);figure(1)%第一幅图plot(x23,y23(:,1),'*r',x45,y45(:,1),'ob',x113,y113(:,1),'+g')%作出各种函数所得结果legend('ode23解','ode45解','ode113解')title('ODE函数求解结果')figure(2)plot(x45,y45)%以ode45为例作出函数以及其各阶导数图legend('y','y一阶导数','y两阶导数')title('y,y一阶导数,y二阶导数函数图')建立求解函数文件myfun03functiondy=myfun03(x,y)dy=zeros(3,1)%初始化变量dy,该行可以没有dy(1)=y(2);%dy(1)表示y的一阶导数,其等于y的第二列值dy(2)=y(3);%dy(2)表示y的二阶导数dy(3)=2*y(3)/x^3+3*y(2)/x^3+3*exp(2*x)/x^3%dy(3)表示y的三阶导数求解过程:[x45,y45]=ode45('myfun03',[1,10],[11030]);第八讲非线性方程求根引言二分法迭代法Newton迭代法MATLAB的非线性方程求根函数小结Matlab中zeroin的算法实现是fzero.
x=fzero(FUN,x0)%x0可以是数,或区间x=fzero(FUN,x0,options)[x,fval]=fzero(FUN,x0,options)[x,fval,exitflag]=fzero(FUN,x0,options)5、MATLAB的非线性方程求根函数例:求方程在x=0.5附近的根。>>x=fzero('x^2+x-1',0.5)x=0.6180例:采用牛顿下山法求方程在区间[1.2,2]上的一个根。>>x=fzero('sqrt(x)-x^3+2',[02])x=1.4759>>r=NewtonDown('sqrt(x)-x^3+2',1.2,2)r=1.4759第九讲解线性方程组的直接解法引言Gauss消元法列主元素消元法误差分析MATLAB的线性方程组求解函数1小结2、Gauss消元法2.1基本思想:逐步消去未知元,将方程组化为与其等价的上三角方程组求解。2.2分两步:第一步:消元过程,将方程组消元化为等价的上三角形方程组;第二步:回代过程,解上三角形方程组,得原方程组的解。2.3Gauss消元的目的a11x1+a12x2+····+a1nxn=b1a21x1+a22x2+····+a2nxn=b2·········································an1x1+an2x2+····+annxn=bn原始方程组约化方程组2.4.1消元过程(化一般方程组为上三角方程组)以四阶为例:其系数增广矩阵为:第一轮消元:计算3个数:[m21
m31
m41]T=[a21
a31
a41]T/a11
用-m21乘矩阵第一行后加到矩阵第二行;
用-m31乘矩阵第一行后加到矩阵第三行;
用-m41乘矩阵第一行后加到矩阵第四行;其系数增广矩阵变为:第二轮消元:计算2个数:[m32
m42]T=[a32(1)
a42(1)]T
/a22(1)
用-m32乘矩阵第二行后加到矩阵第三行;
用-m42乘矩阵第二行后加到矩阵第四行;其系数增广矩阵变为:第三轮消元:计算:m43=a43(2)/a33(2)用-m43乘矩阵第三行后加到矩阵第四行;其系数增广矩阵变为:其对应的上三角方程组为若对于一般的线性方程组Ax=b,其消元过程的计算公式为:(k=1,2,…,n-1)2.4.2回代过程(解上三角方程组)上三角方程组的一般形式为:其中a11…ann≠0回代过程的计算公式:定理:
约化的主元素ak+1,k+1(k)≠0(k=0,1,···,n-1)的充分必要条件是矩阵A的各阶顺序主子式不为零。即注:对角线上的元素ak+1,k+1(k)在Gauss消去法中作用突出,称约化的主元素。推论:
如果A的顺序主子式Dk
≠0(k=1,···,n-1),则Gauss消元法中的约化主元可以表示为例
用高斯消元法求解方程组x1=-13,x2=8,
x3=2m21=3/2m31=4/2m32=-3/0.5>>clearall;A=[234;352;4330];b=[6532]';[L,U]=lu(A);x=U\(L\b)x=-13.00008.00002.0000MATLAB处理:线性代数方程组可以用矩阵形式表示为即:则:5、MATLAB的线性方程组求解函数1a=[0.40960.12340.36780.2943;0.22460.38720.40150.1129;0.36450.19200.37810.0643;0.17840.40020.27860.3927];b=[0.40430.15500.4240-0.2557]’x=a\b或x=v(a)*b第十讲解线性方程组的迭代解法引言简单迭代法赛得尔迭代法迭代解法的收敛性
MATLAB的线性方程组求解函数2小结格式solve('eqn1','eqn2',...,'eqnN','var1,var2,...,varN')5、MATLAB的线性方程组求解函数2格式X=fsolve(FUN,X0)Matlab非线性方程组求解说明:求解方程形式F(X)=0X、F可以是向量或矩阵
X0初值functionF=myfun(x)F=[2*x(1)-x(2)-exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];x=fsolve('myfun',[-1,1])第十一章概率统计函数应用MATLAB提供了大量进行概率统计计算的函数工具,如随机变量的数字特征、特殊分布的概率计算、参数估计函数、假设检验函数、方差分析、回归分析和统计图绘制等。此处主要介绍方差分析和统计图绘制的函数功能。3.1方差分析3.1.1单因素方差分析单因素方差分析是比较两组或多组数据的均值,它返回原假设——均值相等的概率函数anova1格式p=anova1(X)%X的各列为彼此独立的样本观察值,其元素个数相同,p为各列均值相等的概率值,若p值接近于0,则原假设受到怀疑,说明至少有一列均值与其余列均值有明显不同。p=anova1(X,group)%X和group为向量且group要与X对应p=anova1(X,group,'displayopt')%displayopt=on/off表示显示与隐藏方差分析表图和盒图[p,table]=anova1(…)%table为方差分析表[p,table,stats]=anova1(…)%stats为分析结果的构造说明anova1函数产生两个图:标准的方差分析表图和盒图。方差分析表中有6列:第1列(source)显示:X中数据可变性的来源;第2列(SS)显示:用于每一列的平方和;第3列(df)显示:与每一种可变性来源有关的自由度;第4列(MS)显示:是SS/df的比值;第5列(F)显示:F统计量数值,它是MS的比率;第6列显示:从F累积分布中得到的概率,当F增加时,p值
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 【大学课件】数学建模与数学实验 最短路问题
- 2025版农业科技项目种子采购合同范本4篇
- 二零二五年度铝合金户外家具设计与制造合同3篇
- 2025年度新能源设备承包装卸维护合同4篇
- 2025年度橱柜工程承包及质量保证合同4篇
- 二零二五年度足疗养生承包经营合同细则4篇
- 2025年度绿豆农产品绿色生产技术引进合同4篇
- 2025年度塔吊设备承包与绿色施工技术应用合同2篇
- 《物流装卸搬运培训》课件
- 二零二五年度养老产业代售合同协议书4篇
- 2025年度土地经营权流转合同补充条款范本
- 南通市2025届高三第一次调研测试(一模)地理试卷(含答案 )
- 2025年上海市闵行区中考数学一模试卷
- 2025中国人民保险集团校园招聘高频重点提升(共500题)附带答案详解
- 0的认识和加、减法(说课稿)-2024-2025学年一年级上册数学人教版(2024)001
- 重症患者家属沟通管理制度
- 医院安全生产治本攻坚三年行动实施方案
- 法规解读丨2024新版《突发事件应对法》及其应用案例
- 工程项目合作备忘录范本
- 信息安全意识培训课件
- Python试题库(附参考答案)
评论
0/150
提交评论