《数值计算与MATLAB教材》由会员分享，可在线阅读，更多相关《数值计算与MATLAB教材（51页珍藏版）》请在金锄头文库上搜索。
1、数值计算与MATLAB 0 前言 实际问题数学模型数值计算方法程序设计上机计 算求出结果 “数值计算”、“数值分析”为同义词，这门学科有如下特点： u面向计算机：可实现性 u有可靠的理论基础：理论明确 u要有好的计算复杂性：可行性 u可实现数值验证：可验证性 u算法误差明确、并可满足要求：可信性 主要内容：MATLAB基础（方程组、最小二乘法、回归分 析）、MATLAB数据可视化、常用函数、复数、向量计 算、插值法、曲线拟和、数值微积分、多项式、级数展 开、常微分方程求解、符号计算等，其中有些问题包含 在各章节中，主要以MATLAB为线索，进行展开。 1 MATLAB基础 1.1 基本数据类型
2、与运算 MATLAB数据类型及常数表达 名称：字母开头,大小写敏感，缺省变量：ans 实数：科学计数法，特殊数表达eps :浮点数的最 小单位；Inf无穷大；NaN不是一个数值；pi; realmax; realmin 复数：a+bi,或a+bj ，a,b为实数,i、j虚数单位 矩阵：,其中各项之间用空格、逗号、分号 分开，分号换行 1.1 基本数据类型与运算 特殊矩阵：eye ones rand randn zeros (n,m) / (n) :n by m,n by n 线性向量：线性 a:k:b 产生向量 示例：1）用360个弧度表示一圆周； 2）产生200*200个元素为（0，1）正太
3、分布随 机数矩阵； 1.2 运算 参与运算的基本单位：矩阵 单个实数、复数、向量都是矩阵，都可以参加运算 。 矩阵运算： + - * / 数组运算 ： .* ./ . 其中矩阵的除法 / : B/AB*A-1 : ABA-1B 如果XA=B 则:X=B/A 如果 AX=B 则:X=AB 转置 . 非共轭转置 矩阵与常数运算 示例 生成200*200个元素为复数元素矩阵，复数实部与虚部 均为（10，5）正态分布随机数 产生2000个角度为一周均匀分布，直径为正态分布的点 云，方差为10，圆心为（0，0）点。exp 求解 1.3 矩阵基本操作 矩阵的转置 矩阵元素访问及矩阵重构 A（向量，向量）
4、几个常用函数 fliplr flipud rot90 repmat(A,m,n) reshape(A,m,n) 示例： A:100*100, B:40*40,将B嵌入到A中心位置 A:300*400,B:300*400,合成“块交错”的C，块 的大小为5*5，奇数块来至A,偶数块来至于B 1.4 线性方程组、最小二乘法与回归分析 线性方程一般形式 AX=B X=AB A、B可以是复数？ 复数方程，可以做！ A列数与X行数，B 行数一致，未知数 与方程数 A行数少于列数 A行数多于列数 线性方程一般形式AX=B 1）B列数大于1，相当于一个A 对应多个B的方程 2）A行数小于列数，相当于方 程数
5、小于未知数个数 3）A行数大于列数，相当于方 程数多于未知数个数 1.4 线性方程组、最小二乘法与回归分析 未知数多，方程数少：得到一个基本解 未知数少，方程数多：采用最小二乘法，得到一个最综合的 解，也就是使得“=”左端计算的值，与B的差最小。 最小二乘法： 对于 y=ax+b 形式的线性方程，如果已知点(x1,y1)和 (x2,y2)满足方程，则同时满足y1=a*x1+b;y2=a*x2+b，就 可以建立方程 即可求出即可求出( (a,ba,b) )值值, ,如果多于如果多于2 2个点，则应该找到一条直线与点最近，实际上个点，则应该找到一条直线与点最近，实际上 是求出线到点的是求出线到点的
6、y y距离的总和最小，数学上利用距离的总和最小，数学上利用y y距离的平方和确定总误差距离的平方和确定总误差 。即。即E=E=SUM(ySUM(y*-y)2 ,E*-y)2 ,E中仅包含两个未知数中仅包含两个未知数a,ba,b. .求求E E最小可以利用偏导等于最小可以利用偏导等于 0 0的方程。的方程。 1.4 线性方程组、最小二乘法与回归分析 X=1:5 X=x oese(size(x) Y=10 30 32 45 48 A=xy，即为线性方程的两个系数 误差多大？怎么求？（相对误差、 绝对误差） 如何做非线性回归 x=1:10; y=20 21 24 24.6 26 30.1 29.7
7、38 43 53 方程形式为y=ax2+bx+c 实际上 1.4 线性方程组、最小二乘法与回归分析 MATLAB实现: x,y 转换为列向量 X=x.2 x ones(size(x) A=xy, a(1) - a a(2) - b a(3) - c 如何理解ones(size(x)? Size返回什么 其他方程形式回归！ 1.5 环境管理 Clc：清楚控制台（操作窗口） Clear：清楚内存变量，其他形式？清楚一个（一批）特定变 量 help 内容 Quit SAVE fname X Y Z -ASCII -DOUBLe load (filename) who whos 输出格式设定 form
8、at short / long /short e / long e /short g/ long g/ format hex / bank /rat 2 常用函数与数据可视化 2.1 常用函数 Magic(n)魔方阵 pascal(n)阵 hadamard(n) n为2的 密正交阵 calendar(yyyy,mm) 月历 ； clock 时间向量 log2 log10 log sqrt exp abs imag real angle conj fix 向0取整 floor ceil round舍入 sign 数单位值 Sin cos tan asin acos atan atan2 2.2
9、平面曲线绘制 plot(y) 如果y为实数向量则绘制各点折线，复数向量 以实部为x,虚部为y绘制曲线，y为矩阵则以列 为线绘制折线 plot(X1,Y1,.) plot(X1,Y1,LineSpec,.) 绘制（x1,y1）线并指 定线的特性 线型 ：- 实线 -划线 :点线 -.点划线 数据标志：o x + s d * . v x(k+1) a=a+1; else if x(k)0); b=sum(dx5 subp(p+l,0.71*l*exp(i*pi/4); subp(p+l,0.71*l*exp(- i*pi/4); End; end 8 微分方程数值解与数值积分 8.1 微分方程数值
10、 解 可以计算出一个序列的y，即为数值解，可以使用 Cumsum(x)函数实现累加过程，但先要给定初值 。 例如：y=cos(x),y=? - Euler R-K方法微分方程数值解 实际上是前后两点斜率的平均值。上式叫二阶RungeKutta法 上式叫四阶RungeKutta法 MATLAB微分方程数值解法 t，x=ode23(f，ts，x0，options) t，x=ode45(f，ts，x0，options) 分别为二阶，四阶R-K方法，其中f微分方程 （组）的描述函数；ts:为微分点的时间序列 ，求出的x,为与ts对应的解序列，x0为初值 ，options为求解过程误差控制参数，可以省
11、略 示例1 1）定义方程函数 function dx=mf(t,x); dx=cos(t); end 2）求解 ts=0:0.01:10; t,x=ode45(mf,ts,0); Plot(x);%显示曲线 1）定义方程函数 function dy=mf(t,y); dy=-t*y.2; end 2）求解 ts=0:0.01:1; t,y=ode45(mf,ts,2); Plot(y);%显示曲线 示例2 其中：a=20;b=40;c=15; t的范围0,0.5,x(0)=0;y(0)=0 注意：为二元微分方程 1）编辑函数 function dx=mf(t,x) a=20; b=40; c=
12、15; s=sqrt(c-x(1)2+(a*t-x(2)2); dx=b*(c-x(1)/s;b*(a*t-x(2)/s; end 函数返回值dx/dt;dy/dt 2）ts=0:0.01:0.5; x0=0 0; t,x=ode45(mf,ts,x0); %采用缺省的精度求解 X(1)x,x(2)y 示例3 1)定义方程函数 function du=mf(t,u); du=-10*u(1)+9*u(2);10*u(1)-11*u(2); end 2)求解 ts=0:0.01:0.5; t,u=ode45(mf,ts,1 2); 示例4 定义微分方程函数 function dy=mf(t,y)
13、; dy=0;0; dy(1)=y(2); dy(2)=(7*y(2)-5*y(1)+10)/2; end 求解：t,y=ode45(mf,0:0.01:1,2 1) 实际上y(2)，应该与y(1)相同 微分simulink解法 上个方程 8.2 数值积分 trapz(x,y)：梯形法求 积分，求面积 q = quad(fun,a,b) ： Simpson法积分 示例： 定义被积函数 function yy=mf(x); yy=1./(x.3-2*x-5); end 积分运算 Q=quad(mf,0,2) 9 函数极值与优化问题 fplot(fun,limits):绘制函数曲线 如：fplot
14、(sin,0,2*pi); x = fminbnd(fun,x1,x2):返回单变量函数fun， 在x1,x2区间上的 最小值， 如：fminbed(sin,-1,1),结果为-1 x = fminsearch(fun,x0)：返回多变量函数fun，在x0附近的最小值 如：function f=myf(x); f=100*(x(2)-x(1).2).2+(1-x(1).2; end x=fminsearch(myf,0 1) 9 函数极值与优化问题 x = bintprog(f, A, b, Aeq, beq)：求多变量函数线性函数f， (1/0)规划的极小值，约束条件：AX=b;AeqX=b
15、eq;X取值 范围1/0,其中f为向量表示的线性目标函数 如：求f(x)=-9x1-5x2-6x3-4x4极小值，满足 6x1+3x2+5x3+2x4=9; x3+x4=1; * x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)：目 标规划 x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub)：条件极值 9函数极值与优化问题 x = fsolve(fun,x0, options) ：返回fun在x0的解，如： 在-5 -5附近的解，定义函数 function f=myf1(x) f=2*x(1)-x(2)-exp(-x(1); -x(1)+2*x(2)-exp(-x(2); end 求解：x=fsolve(myf1,-5;-5) options = optimset(param1,value1,param2,value2,.) Display:Level of display. off displays no output; iter displays output at e
