Matlab初步（讲稿）
上课方式：学生边听讲、边用机器练习。
调用Matlab软件：在windows平台上，双击“Matlab”图标。
说明：执行此命令，将进入“Matlab工作区（命令区）”，在这里可以下达、执行符合Matlab语法的各种命令。
一．矩阵；数组
1．矩阵
例：输入一个矩阵，并用A代表此矩阵；
再输入一个矩阵，并用a代表此矩阵,
命令为：A=[2,-18;0,31;-59,6]; a=[0,1,-1;-1,2,1];
说明：（1）输入上述命令后，敲回车键，机器才执行此命令
（2） 你发现机器屏幕没反应，其实它早已接受并执行了你的命令，不信？ 请下命令 A 就有反应了。 再下命令 a
（3）清屏（把屏幕上的Matlab工作区清理的干干净净）
屏幕空了，但你的那两个矩阵A、a还在机器中。
（4）注意要点：中括号 逗号 分号 字母大小写。
注意：（1） A(i,j) 表示矩阵A的第i行、第j列交叉处的元素练习,A A(3,1) A(1,2) A(2,3) a a(2,3) y=A(3,1)*a(2,2)
（2）可修改个别元素 练习,a(2,2)=8 a
（3）可把矩阵的型号（即：行个数、列个数）放大
练习,A(2,4)=9 A a(3,5)=99 a
（4）一些特殊矩阵
m行n列的 全0矩阵：zeros(m,n) 全1矩阵：ones(m,n)
m行n列的 单位矩阵：eye(m,n) 随机矩阵：rand(m,n)
（随机矩阵的每个元素都是：开区间(0.1)内的均匀分布随机数）
练习,zeros(2,3) zeros(1,5) ones(4,4) ones(2,4)
eye(3,3) eye(3,5) eye(4,2) rand(3,4) rand(1,10)
请产生20个在区间(0,8)内的随机数,请产生20个在区间(3,8)内的随机数
2．对矩阵作裁剪、拼接裁剪：从矩阵中提取某些行、某些列 （关键符号,）
如（练习）,A(2,:) 是A的第2行 A(:,1) 是A的第1列
A(1:2,2:4) 是A的第1、2行，与第2、3、4列交叉点元素
输入 ，它的第2、3行，与第3、4、5、6列交叉点元素是什么?
W的第1、3、5、7列构成的矩阵是 W(:,1:2:7) （起点1，步长2，终点7）
W的第1、3行，第2、5、8列构成的矩阵是 W(1:2:3,2:3:8)
问：W的第1、4、7列构成的矩阵？ W的第1、3行，第2、4、6列构成的矩阵？
检验,W(:,7:-2:1) 看结果，猜一猜什么规则？
检验,W W(:,5)=[] 看结果，猜一猜什么规则？
此时，W还剩7个列，请你用一个命令去掉它的第3、5列。
拼接：把若干个矩阵、数，拼凑、结合成一个矩阵先做准备：把机器中全部的常量（包括矩阵、数）清除，命令为 clear
再重新输入：，
练习,左右拼接 [A,B] [B,A] [a,8] [8,a,7,6]
上下拼接 [A;a] [a;A] [A;8,18,28;a]
试一试,左右拼接 [A,a] [A,8]
上下拼接 [A;B] [A;8,18]
题：用这三个矩阵A、B、a，拼接出矩阵
3．数组
型矩阵（即：只有一行的矩阵），称为行向量，也称为数组。
例：用a1表示一个从0到18的全体偶数构成的数组。
解,法一,格式为 起点,步长,终点
命令为 a1=0:2:18
法二,格式为 linspace(起点,终点,个数)
命令为 a1=linspace(0,18,10)
（这两种格式的命令都表示等差数列）
练习：用a2表示一个从9到1的全体奇数构成的数组；
用a3表示一个从9到-9的全体整数构成的数组；
用第二种格式、用a4表示一个从-6到8的5个数构成的等差数列；
把两个数组a2与a4合并成一个数组； [a2,a4]
用a5表示数组 1,4,7,…,37,38,35,32,…,2,3,6,9,…,39
二．运算
1．矩阵运算加 减 乘 乘方 左除 右除 转置
+ - * ^ \ / ‘ (单引号)
注：（1） +,-,*,^,\,/ 都应符合矩阵运算规则；
（2）特别，数与数的运算,加 + 减 - 乘 * 除 /
（3）数与矩阵的 加 减 乘练习,A 2+A 2-A A-2 2*A A*2
（4） A\B 读作 A左除B,其本质是
A/B 读作 B右除A,其本质是
主要用来解“矩阵方程”，如：AXB=C，其中A，B，C 均为已知矩阵，X是未知矩阵，则 ,命令是 X=(A\C)/B,
例：，求X=？
例：解线性方程组
2．数组运算点乘 点乘方 点左除 点右除
,*,^,\,/
练习,a=[1,2,3,4] b=[5,6,7,8] c=[5,6,7]
a.*b b.^2 a.^3 a.^0.5 a.\b a./b b.\a a.*c
通过练习，搞明白运算规则。
若两个矩阵A、B同型号，则 A.*B A.^B 合法；
若两个矩阵A、B同型号，且B的元素都非零，则 A./B B.\A 合法,
练习：，学习运算规则。
三．命令语句，函数
1．命令语句格式,变量=表达式
（1）“表达式”不可缺省,“变量=”可缺省。
若语句中缺省“变量=”，则机器把执行结果自动记为 ans,
例如：[1,0,-2]*[3;1;2],(此命令中只有表达式，缺省“变量=”)，结果为 ans = -1,
（2）同一行可以写多个语句，语句与语句之间用逗号或分号隔开。
练习：（借用前面的a,b）
一整行命令为 a,b,c=a*b’,d=a.*b,e=a./b 再敲回车键执行命令。
（3）若只要求机器执行命令、而不让机器把某个语句的执行结果显示在屏幕上,则必须在该语句之后跟分号。
练习：一整行命令为 a,b;c=a*b’,d=a.*b;e=a./b; 再敲回车键执行命令。
（4）关于变量名：你可以任意用一个字符串来代表一个变量名，但必须满足两条：字母打头；避开Matlab软件的专用符号。
几个专用符号,pi 是圆周率
eps 是最小浮点数（即：机器所能够处理的最小正数）
Inf 是
NaN 是不定值练习,pi eps 2/0 (3-3)/(2*0)
介绍 format long 与 format
2．函数
（1）普通函数（也称标量函数，简称函数）
常见的函数,sin cos tan exp log log10 sqrt abs
round floor ceil
四舍五入取整 负向取整 正向取整设是普通函数，，则
练习,，cos(A) tan(A) exp(A) 1og(A) log10(A) sqrt(A)
a=[-6.01 -4.49 -2.50 -0.99 0.99 2.50 4.49 6.01]
round(a) floor(a) ceil(a)
（2）向量函数（也称数组函数）
常见的向量函数,
length max min sum sort prod mean median
长度(维数) 和 排序 乘积 平均值 中值定义：设a是一个行向量、或列向量，则
length(a) 是指 向量a所包含的数的个数；
max(a) 是指 向量a中的最大数；
sum(a) 是指 向量a中的全部数的和；
sort(a) 是把向量a中的数，按照从小到排序，得到一个新的向量；
prod(a) 是指 向量a中的全部数的积；
mean(a) 是指 向量a中的全部数的平均值；
median(a) 若length(a)奇，则是指 向量sort(a)中位于中间的那个数；
若length(a)偶，则是指 向量sort(a)中位于中间的那两个数的平均值。
例：,则 sort(a)=[-1,0,3,4,6],sort(b)=,
median(a)=3,median(b)=0.5,mean(a)=2.4,mean(b)=2.25
继续定义：设A是一个真正的矩阵（行数、列数都不小于2），则
length(A) 是A的行数、列数两者中的最大者；
sort(A) 是 把A的每个列各自排序，得到一个新的矩阵；
其余的6个向量函数，作用于A 是 首先分别作用于每个列得到一个个数，
然后再把这些数写成一行，得到一个行向量。
定义完毕。
例：，则 length(A)=5 max(A)=[6,8,5,5]
min(A)=[0,-5,-4,1] sum(A)=[12,3,4,15] sort(A)=
mean(A)=[2.4,0.6,0.8,3] median(A)=[2,0,1,3]
习题：用机器产生一个型矩阵G，它的元素是“0至99范围内的均匀分布随机自然数”，G的第3行记为g1，G的第4列记为g2，用机器计算G,g1,g2 的所有八种向量函数值。
（3）矩阵函数第一类矩阵函数：（用来产生特殊矩阵）
如：zeros ones eye rand randn（标准正态分布：N(0.1)）等等。
第二类矩阵函数：（用来计算矩阵的 行列式、逆、秩、特征值 等等）
如,det inv rank eig poly 等等。
行列式 逆 秩 特征值 特征多项式例：B=,计算B的行列式、逆、秩、特征值、特征多项式。
（课堂上教师详细做。讲解poly：多项式系数（从高次到低次）构成的向量）
四．多项式例：用Matlab表示多项式
解,P=[1,-6,0,2,-8];
Px=poly2str(P,’x’)
执行结果为 Px=x^4-6*x^3+2*x-8
说明：（1）一个多项式的本质核心是它的各次项系数，把这些系数写成一个行向量（从高次到低次）称为该多项式的系数向量，系数向量完全决定了多项式，上面第一个命令就是把系数向量记作变量P；
（2）poly2str(P,’x’) 就是让机器根据系数向量P、以及符号’x’写出多项式。若把此句改为PPP=poly2str(P,’t’)，则执行结果为PPP=t^4-6*t^3+2*t-8 。
习题：用Matlab表示多项式
例：已知6次多项式G(x)的6个根为
-0.38 0.5 4 4.77 0.3+0.4i 0.3-0.4i
请写出G(x)的表达式。
解：R=[-0.38,0.5,4,4.77,0.3+0.4i,0.3-0.4i];
P=poly(R); （根据根R计算出系数向量）
P1=real(P); （因机器计算有误差，P带有虚部。根据数学知识，当虚根共轭成对出现时，必为实系数多项式。指令real(P)就是提取P的实部）
Gx=polt2str(P1,’x’)
习题：5次多项式F(x)的5个根为 -2 -1 0 1 2,请写出G(x)的表达式。
7次多项式H(x)的7个根为 3.2,-1.8,1,i,32i,请写出H(x)的表达式。
小结：指令poly(A)，若A为方阵，则产生A的特征多项式的系数；
若A为行向量，则产生以A为根的多项式的系数。
例：
解,命令,
再命令：c=conv(a,b); (产生系数向量)
Cx=polt2str(c,’x’)
完毕。
习题：=？
指令 roots(P) 的功能是：计算出以P为系数向量的多项式的根。
例：（求的根）
指令 polyval(P,A) 的功能是：计算出以P为系数向量的多项式在A的每个元素处的值。
例：
,求
解：命令为 P=[6,-8,0,3,-2,1];B=[2,-4,6;-3,5,-7];polyval(P,B) 完毕。
指令 polyvalm(P,A) 的功能是：A是方阵，按照矩阵运算规则计算出矩阵多项式的值。
例：
,求
解：命令为 P=[6,-8,0,3,-2,1];A=[2,-4;-3,5];polyvalm(P,A) 完毕。
五．Matlab工作区下的操作（初级）
清屏,clc
显示全部变量名,who（简显） 或 whos（详细显示）
显示某个变量AAA的内容,AAA 或 disp(AAA)
清除变量：（1）只清除变量AAA，则命令是 clear(AAA)
（2）清除机器中的全部变量，则命令是 clear
数据的输出格式,format 短型
format long 长型
format rat 分数型六．图形功能
1．二维图形 （命令词,plot ）
例1：已知x取值为 1,2,3,4,5,6
y取值为 0,0.5,0.7,1.1,0.9,0.2
请以x为横坐标、y为纵坐标，画出折线图。
解：以下三个方法中的任意一个都可以，其结果完全一样
（1）x=[1,2,3,4,5,6];y=[0,0.5,0.7,1.1,0.9,0.2];plot(x,y)
（2）x=1:6;y=[0,0.5,0.7,1.1,0.9,0.2];plot(x,y)
（3）y=[0,0.5,0.7,1.1,0.9,0.2];plot(y)
例2：画出正弦曲线y=sin(x)在区间上的图形。
练习： 画图一图多线：（在同一个坐标图中画两条以上的曲线）
法一：例1：x=0:pi/10:2*pi;y1=sin(x);y2=cos(x);plot(x,y1,x,y2)
例2：x=0:pi/10:2*pi;y1=sin(x);y2=cos(x);
y3=(sin(x)+cos(x))/2;plot(x,y1,x,y2,x,y3)
例3：在同一图中画四条线，横坐标范围是[0，2]，函数分别是
y= y= y= y
x1=0:pi/10:2*pi;y1=(1+sin(x))/2;x2=[0,2*pi];y2=x2/(2*pi);
y3=x1.*(2*pi-x1)/pi^2;y4=[0.6,0.6];
plot(x1,y1,x2,y2,x1,y3,x2,y4)
法二：重做例3：只需把plot(x1,y1,x2,y2,x1,y3,x2,y4)改成
hold on,plot(x1,y1),plot(x1,y3),plot(x2,y2),plot(x2,y4),hold off
习题：在同一图中画四条线，横坐标范围是[1，2+1]，函数分别是
,,,
加“坐标网格”：
法一：在 plot(……) 之后紧跟 grid
演示：plot(x1,y1,x2,y2,x1,y3,x2,y4),grid
或者hold on,plot(x1,y1),plot(x1,y3),plot(x2,y2),plot(x2,y4),grid,hold off
法二：事后补加（先画好没有网格的图，再补加网格）
演示：先命令plot(x1,y1),再命令grid
加“标记”：（1）给图加标题：title(‘******’)
（2）给坐标轴加标记：横轴：xlabel(‘******’)
纵轴：ylabel(‘******’)
（3）给曲线加标记：text(aaa,bbb,‘******’)
说明：这四个命令的用法与前面grid相同（紧跟plot(….)，或事后补加）；
红色符号是本质（包括：小括号、逗号、单引号），不可改变；
***** 是你准备加注的字符，由你自己决定；
你准备把“曲线标记”放在哪个位置？ 答：aaa是横坐标，bbb是纵坐标。
演示：略写，教师在课堂机器上演示。
对坐标系进行控制：
axis equal 使得横轴与纵轴的单位长度相同；
axis square 使得图呈现正方形；
axis(aax,bbx,aay,bby) 实现你的愿望：使得横坐标从aax到bbx、纵坐标从aay到bby,
演示：略写，教师在课堂机器上演示。
多幅图形：在同一个画面上建立多个坐标系，每个坐标系中各自画图、互不干扰关键词：subplot(m,n,p),把画面分割成个图形区域，p代表第几个图。
例如：借用前面的 x1,x2,y1,y2,y3,y4,
subplot(1,4,1),plot(x1,y1)
subplot(1,4,2),plot(x2,y2),grid
subplot(1,4,3),plot(x1,y3),grid
subplot(1,4,4),plot(x2,y4)
结果真难看！ 改
subplot(2,2,1),plot(x1,y1)
subplot(2,2,2),plot(x2,y2),grid
subplot(2,2,3),plot(x1,y3),grid
subplot(2,2,4),plot(x2,y4)
结果好看了！
习题：同一个画面上6幅图：
2．三维图形
先不动脑筋，看别人做的一个三维曲面图，感受一下Matlab作图效果。
例：画出曲面
解：x=-7.5:0.5:7.5;y=x; [X,Y]=meshgrid(x,y);
R=sqrt(X.^2+Y.^2)+eps;Z=sin(R)./R;
mesh(X,Y,Z)
（1）带网格的曲面图 (关键词：mesh )
例如：画曲面
x=-0.5:0.1:0.5;y=x;
[X,Y]=meshgrid(x,y); （组建三维图形的X,Y数组）
Z=sqrt(1-X.^2-Y.^2);
mesh(X,Y,Z) （画图）
（2）曲面的等高线图 (关键词：contour3 )
续上题：把 mesh(X,Y,Z) 改为 contour3(X,Y,Z,5)
(其中，5 是准备画5条等高线)
（3）空间三维曲线图 (关键词：plot3 )
例：螺旋线的参数方程为：
试在范围内，画出该曲线。
解：t=-20:pi/30:20;x=sin(t);y=cos(t);z=t;plot3(x,y,z)
习题：曲面，分别画网格图和等高线图。
曲线 自己适当取t的范围。
七．程序
1．编写、保存、打开、修改、执行程序编写新程序：击活菜单File，选New，单击M—file，进入编辑区，可编写程序。
例1：（演示：怎样编写、保存、打开、修改、执行程序？）
请按照规律1+2=3，2+3=5，3+5=8，5+8=13写出30个数作为一个数组a=[1,2,3,5,8,13,……],
简要说明：（1）保存新文件：File----Save,文件名必须字母打头；
（2）打开旧文件：File----Open----文件名，打开以后就可修改内容；
（3）执行：Debug----Run,
2．关系运算，逻辑运算
（1）关系运算：（比较两个数）
共6种运算符,< > <= >= == ~=
小于 大于 小于或等于 大于或等于 等于 不等于
练习,1<2 1>2 2<2 2<=2 1==2 1~=2
课堂判断,？？？？？？
Matlab规定：当关系成立时，结果是1；当关系不成立时，结果是0,
课堂做，且机器检验,(1<2)=? (1>2)=? (2<2)=?
(2<=2)=? (1==2)=? (1~=2)=?
习题,a=1:5 b=5:-1:1
c=a>=3,d=a==3,e=a<b,f=a==b
先手工算出 数组c,d,e,f,再上机器检验。
（2）逻辑运算,共有3种,& | ~
与(并且) 或 非(不是)
Matlab规定：(1) 当关系真时，结果是1；当关系假时，结果是0 ；
（2）若关系运算和逻辑运算混合了，则，关系运算优先，逻辑运算靠后，整体还叫关系表达式。
学习体会：这个符号“=”，它不叫“等号”，而叫“赋值号”。关系符号、逻辑符号都比它优先。
练习：a=1:8 b=8-a c=~(a>4) d=(a>=3)&(b<6)
则 a=
b=
c=
d=
3．选择结构（又称：条件语句）
格式一,if <关系表达式> 若 关系表达式 的值是1，则执行语句1；
<语句1> 若 关系表达式 的值是0，则不执行语句1，而去
end 执行end后面的语句。
格式二,if <关系表达式> 若 关系表达式的值是1，则执行语句1,再跳到end后面；
<语句1> 若 关系表达式的值是0，则不执行语句1，而执行语句2
else
<语句2>
end
格式三,if <关系表达式1> 若 关系表达式1的值是1，
<语句1> 则执行语句1,再跳到end后面；
elseif <关系表达式2> 若 式1假且式2真，
<语句2> 则执行语句2,再跳到end后面；
elseif <关系表达式3> 若 式1假、式2假且式3真，
<语句3> 则执行语句3,再跳到end后面；
…… …… …… 略。
end
4．循环结构格式一,for <循环参数> = <初值>,<步长>,<终值>
<语句1>
end
说明：对于每个参数值，都重复执行一次语句1,(课堂上，借用前面程序作解释)
例2：请输入、并显示一个20阶方阵，它的元素满足
解,笨法一：A=[1,2,3,4,…,19,20;2,3,4,…,21;……]
方法二：编程，程序清单如下：
for i=1:20
for j=1:20
A(i,j)=i+j-1;
end
end
A
( 循环结构中套循环！)
有时，程序中必需循环结构，但找不道合适的循环参数，或不知道循环次数，就用下面格式：
格式二,while <关系表达式1>
<语句1>
end
说明：每次执行完语句1，都要倒回去检验关系表达式1，只要关系表达式1的值是1，就再执行语句1。就这样重复执行直到关系表达式1的值是0为止。初学者慎用此结构，避免陷入死循环。
死循环例子：a=0;
while a>=0
a=a+1;
end
例3：记，求自然数m，使得
解,(课堂现编)
5．关于函数M----文件格式,function <被输出的变量名>=<文件名>（全部被引入的变量名）
注：变量名之间用逗号隔开例如，函数文件 jiecheng.m 的程序内容为：
function J=jiecheng (n)
J=1;
for k=1:n
J=J*k;
end
准备好该文件后，回到Matlab命令区。此时，你想计算几的阶乘？ 想算 3!，5!，12!，100!，只需分别下命令 jiecheng (3),jiecheng (5),jiecheng (12),jiecheng (100) 即可。
所以，文件jiecheng.m就是计算n阶乘的函数文件。
八．举例（用Matlab作计算）
小题1：随机产生一个满足“标准正态分布”的100维数组，记作a,再把数组a的这100个数的次序完全颠倒，得到一个新数组，记作b，输出b,
小题2：随机产生一个满足“(0,1)内均匀分布”的型矩阵，求出其最大元素及其所处的位置。
小题3：随机产生一个满足“(1,10)内均匀分布”的20维数组，用起泡法对20个数由小到大排序,即将相邻两个数比较,将小的调到前头.（不准用sort命令）。
小题4：编程求
小题5,一球从100米高度自由落下,每次落地后反跳回原高度的一半,再落下,求它在第10次落地时,共经过多少米?第10次反弹有多高?
大题1：
汽车头部的车灯，其反射面为一个旋转抛物面，方程为
。
灯丝（即：线光源）是一个直线段、长度为4、位于Y轴上且以原点为中点。 设点P(0,t,0)是线光源上任一点，则 设点M()是旋转抛物面任一点。从P到M的射线的反射线与平面 的交点记为N(25000,y,z)，计算公式为
点P在线光源上遍历，点M在抛物面上遍历，则，对应的点N在平面 上移动轨迹产生的区域就是反射光亮区。取t以步长0.2从-2到2，以步长1从-14.4到6.6。每固定一个，分别再遍历以步长从0到,编程画出反射光亮区。
大题2：
x
0
3
5
7
9
11
12
13
14
15
y
0
1.2
1.7
2.0
2.1
2.0
1.8
1.2
1.0
1.6
用Lagrange插值法，求当x分别取0,0.2,0.4,0.6,0.8,……,15时，对应的函数值y，并画出函数y=y(x)的图象.
提示：Lagrange插值公式
,
大题3,已知：，测试得下面数据
x,0 1153 2045 2800 3466 4068 4621 5135 5619 6152
y,0 20 40 60 80 100 120 140 160 183.5,
试用“最小二乘法”估计参数，并算出节点处的总误差.
大题4：（数值积分的龙贝格算法）计算，解法如下：
(k=1,2,3,……)
每当算出一个，就立即按下面公式计算
计算，输出格式为 format long,算到为止。顺便看看与准确值有多大误差？
大题5,求方程的最小正根。
要求：先观察确定有根区间[a,b]，区间宽度不大于0.2 ；
以下工作由程序来做：用“二分法”把有根区间宽度缩小到0.01，再启动“牛顿迭代公式”计算10次，输出近似根，输出格式为 format long
提示：0是整数？（笔误，改正：0是正数？）
方程f(x)=0,一个近似根记为，求根的牛顿迭代公式为
,
九．补充内容：
用Matlab软件求解微分方程
1．解析解
（1）一阶微分方程求的通解：dsolve('Dy=1+y^2','x')
求的通解：dsolve('Dy=1+x^2-y','x')
求的特解：dsolve('Dy=1+y^2',’y(0)=1’,'x')
（2）高阶微分方程求解其中，，命令为：
dsolve('x^2*D2y+x*Dy+(x^2-0.5^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x')
求的通解，命令为：
dsolve('D3y-2*Dy+y-4*x=0','x')
输出为：
ans=8+4*x+C1*exp(x)+C2*exp(-1/2*(5^(1/2)+1)*x)+C3*exp(1/2*(5^(1/2)-1)*x)
（3）一阶微分方程组求的通解：[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','x')
输出为,f =exp(3*x)*(cos(4*x)*C1+sin(4*x)*C2)
g =-exp(3*x)*(sin(4*x)*C1-cos(4*x)*C2)
若再加上初始条件，则求特解：
[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1','x')
输出为,f =exp(3*x)*sin(4*x)
g =exp(3*x)*cos(4*x)
2．数值解
（1）一阶微分方程
现以步长h=0.1用“4阶龙格—库塔公式”求数值解：
先建立“函数M—文件”：function f=eqs1(x,y)
f=y-2*x/y;
再命令,格式为：
[自变量,因变量]=ode45(‘函数文件名’,节点数组,初始值)
命令为,[x,y]=ode45('eqs1',0:0.1:1,1)
若还要画图，就继续命令,plot(x,y)
（2）一阶微分方程组
只须向量化，即可用前面方法：
function f=eqs2(x,y)
f=[cos(x)+2*y(1)-y(2);sin(x)-y(1)+2*y(2)];
将此函数文件，以文件名eqs2保存后，再下命令：
[x,y]=ode45('eqs2',0:0.1:1,[0.2;0.3])
（注：输出的y是矩阵，第i列为函数的数值解）
要画图，继续命令：hold on,plot(x,y(:,1)),plot(x,y(:,2)),hold off
（3）高阶微分方程先化成一阶微分方程组，再用前面方法。
上机练习：
准备：令，化成
用机器,函数文件eqs3内容？
命令？
画图？
用Matlab作最小二乘曲线拟合
1．用n次多项式作最小二乘拟合已知，要从（即：全体次数不高于n的多项式集合）中找一个，使得在节点处的总误差达到最小。
Matlab命令格式：系数数组=polyfit(节点数组，函数值数组，次数n)
例1：对函数C=C(t)测量得下面一组数据：
t,1 2 3 4 5 6 7 8 9
C：4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50
试分别用1次、2次、6次多项式作拟合，并画图显示拟合效果。
clear
hold on
x0=1:9;y0=[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50];
for i=1:9
plot(x0(i),y0(i),'+')
end
a1=polyfit(x0,y0,1),a2=polyfit(x0,y0,2),a6=polyfit(x0,y0,6)
x=0:0.1:10;
y1=polyval(a1,x);y2=polyval(a2,x);y6=polyval(a6,x);
plot(x,y1,x,y2,x,y6)
hold off
为了准确判断拟合效果，需计算“节点处的总误差”：（续前面程序）
wc1=sqrt(sum((polyval(a1,x0)-y0).^2))
wc2=sqrt(sum((polyval(a2,x0)-y0).^2))
wc6=sqrt(sum((polyval(a6,x0)-y0).^2))
2．用一般函数作最小二乘拟合已知，要用一个函数来近似代表，此函数中含有几个待定参数，现在的任务是：确定参数的值，使得在节点处的总误差达到最小。
用Matlab计算：先按照建立一个函数文件，文件第一句格式为：
function f=文件名(参数组,节点数组)
再用下面命令格式：
参数组=lsqcurvefit(‘函数文件名’，参数组初值，节点数组，函数值数组)
例2：对例1中的数据，用函数作拟合，并画图显示拟合效果，给出节点处总误差。
先建立并保存函数文件：文件名syp78hswj 内容为：
function f=syp78hswj(a,x0)
f=a(1)+a(2)*exp(-a(3)*x0);
再下面主程序：
clear
hold on
x0=1:9;y0=[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50];
for i=1:9
plot(x0(i),y0(i),'+')
end
cscz=[1,1,1];
a=lsqcurvefit('syp78hswj',cscz,x0,y0)
x=0:0.2:10;f=syp78hswj(a,x);plot(x,f)
f=syp78hswj(a,x0);
wc=sqrt(sum((f-y0).^2))
hold off
执行得：
节点处总误差 wc=0.0076
上课方式：学生边听讲、边用机器练习。
调用Matlab软件：在windows平台上，双击“Matlab”图标。
说明：执行此命令，将进入“Matlab工作区（命令区）”，在这里可以下达、执行符合Matlab语法的各种命令。
一．矩阵；数组
1．矩阵
例：输入一个矩阵，并用A代表此矩阵；
再输入一个矩阵，并用a代表此矩阵,
命令为：A=[2,-18;0,31;-59,6]; a=[0,1,-1;-1,2,1];
说明：（1）输入上述命令后，敲回车键，机器才执行此命令
（2） 你发现机器屏幕没反应，其实它早已接受并执行了你的命令，不信？ 请下命令 A 就有反应了。 再下命令 a
（3）清屏（把屏幕上的Matlab工作区清理的干干净净）
屏幕空了，但你的那两个矩阵A、a还在机器中。
（4）注意要点：中括号 逗号 分号 字母大小写。
注意：（1） A(i,j) 表示矩阵A的第i行、第j列交叉处的元素练习,A A(3,1) A(1,2) A(2,3) a a(2,3) y=A(3,1)*a(2,2)
（2）可修改个别元素 练习,a(2,2)=8 a
（3）可把矩阵的型号（即：行个数、列个数）放大
练习,A(2,4)=9 A a(3,5)=99 a
（4）一些特殊矩阵
m行n列的 全0矩阵：zeros(m,n) 全1矩阵：ones(m,n)
m行n列的 单位矩阵：eye(m,n) 随机矩阵：rand(m,n)
（随机矩阵的每个元素都是：开区间(0.1)内的均匀分布随机数）
练习,zeros(2,3) zeros(1,5) ones(4,4) ones(2,4)
eye(3,3) eye(3,5) eye(4,2) rand(3,4) rand(1,10)
请产生20个在区间(0,8)内的随机数,请产生20个在区间(3,8)内的随机数
2．对矩阵作裁剪、拼接裁剪：从矩阵中提取某些行、某些列 （关键符号,）
如（练习）,A(2,:) 是A的第2行 A(:,1) 是A的第1列
A(1:2,2:4) 是A的第1、2行，与第2、3、4列交叉点元素
输入 ，它的第2、3行，与第3、4、5、6列交叉点元素是什么?
W的第1、3、5、7列构成的矩阵是 W(:,1:2:7) （起点1，步长2，终点7）
W的第1、3行，第2、5、8列构成的矩阵是 W(1:2:3,2:3:8)
问：W的第1、4、7列构成的矩阵？ W的第1、3行，第2、4、6列构成的矩阵？
检验,W(:,7:-2:1) 看结果，猜一猜什么规则？
检验,W W(:,5)=[] 看结果，猜一猜什么规则？
此时，W还剩7个列，请你用一个命令去掉它的第3、5列。
拼接：把若干个矩阵、数，拼凑、结合成一个矩阵先做准备：把机器中全部的常量（包括矩阵、数）清除，命令为 clear
再重新输入：，
练习,左右拼接 [A,B] [B,A] [a,8] [8,a,7,6]
上下拼接 [A;a] [a;A] [A;8,18,28;a]
试一试,左右拼接 [A,a] [A,8]
上下拼接 [A;B] [A;8,18]
题：用这三个矩阵A、B、a，拼接出矩阵
3．数组
型矩阵（即：只有一行的矩阵），称为行向量，也称为数组。
例：用a1表示一个从0到18的全体偶数构成的数组。
解,法一,格式为 起点,步长,终点
命令为 a1=0:2:18
法二,格式为 linspace(起点,终点,个数)
命令为 a1=linspace(0,18,10)
（这两种格式的命令都表示等差数列）
练习：用a2表示一个从9到1的全体奇数构成的数组；
用a3表示一个从9到-9的全体整数构成的数组；
用第二种格式、用a4表示一个从-6到8的5个数构成的等差数列；
把两个数组a2与a4合并成一个数组； [a2,a4]
用a5表示数组 1,4,7,…,37,38,35,32,…,2,3,6,9,…,39
二．运算
1．矩阵运算加 减 乘 乘方 左除 右除 转置
+ - * ^ \ / ‘ (单引号)
注：（1） +,-,*,^,\,/ 都应符合矩阵运算规则；
（2）特别，数与数的运算,加 + 减 - 乘 * 除 /
（3）数与矩阵的 加 减 乘练习,A 2+A 2-A A-2 2*A A*2
（4） A\B 读作 A左除B,其本质是
A/B 读作 B右除A,其本质是
主要用来解“矩阵方程”，如：AXB=C，其中A，B，C 均为已知矩阵，X是未知矩阵，则 ,命令是 X=(A\C)/B,
例：，求X=？
例：解线性方程组
2．数组运算点乘 点乘方 点左除 点右除
,*,^,\,/
练习,a=[1,2,3,4] b=[5,6,7,8] c=[5,6,7]
a.*b b.^2 a.^3 a.^0.5 a.\b a./b b.\a a.*c
通过练习，搞明白运算规则。
若两个矩阵A、B同型号，则 A.*B A.^B 合法；
若两个矩阵A、B同型号，且B的元素都非零，则 A./B B.\A 合法,
练习：，学习运算规则。
三．命令语句，函数
1．命令语句格式,变量=表达式
（1）“表达式”不可缺省,“变量=”可缺省。
若语句中缺省“变量=”，则机器把执行结果自动记为 ans,
例如：[1,0,-2]*[3;1;2],(此命令中只有表达式，缺省“变量=”)，结果为 ans = -1,
（2）同一行可以写多个语句，语句与语句之间用逗号或分号隔开。
练习：（借用前面的a,b）
一整行命令为 a,b,c=a*b’,d=a.*b,e=a./b 再敲回车键执行命令。
（3）若只要求机器执行命令、而不让机器把某个语句的执行结果显示在屏幕上,则必须在该语句之后跟分号。
练习：一整行命令为 a,b;c=a*b’,d=a.*b;e=a./b; 再敲回车键执行命令。
（4）关于变量名：你可以任意用一个字符串来代表一个变量名，但必须满足两条：字母打头；避开Matlab软件的专用符号。
几个专用符号,pi 是圆周率
eps 是最小浮点数（即：机器所能够处理的最小正数）
Inf 是
NaN 是不定值练习,pi eps 2/0 (3-3)/(2*0)
介绍 format long 与 format
2．函数
（1）普通函数（也称标量函数，简称函数）
常见的函数,sin cos tan exp log log10 sqrt abs
round floor ceil
四舍五入取整 负向取整 正向取整设是普通函数，，则
练习,，cos(A) tan(A) exp(A) 1og(A) log10(A) sqrt(A)
a=[-6.01 -4.49 -2.50 -0.99 0.99 2.50 4.49 6.01]
round(a) floor(a) ceil(a)
（2）向量函数（也称数组函数）
常见的向量函数,
length max min sum sort prod mean median
长度(维数) 和 排序 乘积 平均值 中值定义：设a是一个行向量、或列向量，则
length(a) 是指 向量a所包含的数的个数；
max(a) 是指 向量a中的最大数；
sum(a) 是指 向量a中的全部数的和；
sort(a) 是把向量a中的数，按照从小到排序，得到一个新的向量；
prod(a) 是指 向量a中的全部数的积；
mean(a) 是指 向量a中的全部数的平均值；
median(a) 若length(a)奇，则是指 向量sort(a)中位于中间的那个数；
若length(a)偶，则是指 向量sort(a)中位于中间的那两个数的平均值。
例：,则 sort(a)=[-1,0,3,4,6],sort(b)=,
median(a)=3,median(b)=0.5,mean(a)=2.4,mean(b)=2.25
继续定义：设A是一个真正的矩阵（行数、列数都不小于2），则
length(A) 是A的行数、列数两者中的最大者；
sort(A) 是 把A的每个列各自排序，得到一个新的矩阵；
其余的6个向量函数，作用于A 是 首先分别作用于每个列得到一个个数，
然后再把这些数写成一行，得到一个行向量。
定义完毕。
例：，则 length(A)=5 max(A)=[6,8,5,5]
min(A)=[0,-5,-4,1] sum(A)=[12,3,4,15] sort(A)=
mean(A)=[2.4,0.6,0.8,3] median(A)=[2,0,1,3]
习题：用机器产生一个型矩阵G，它的元素是“0至99范围内的均匀分布随机自然数”，G的第3行记为g1，G的第4列记为g2，用机器计算G,g1,g2 的所有八种向量函数值。
（3）矩阵函数第一类矩阵函数：（用来产生特殊矩阵）
如：zeros ones eye rand randn（标准正态分布：N(0.1)）等等。
第二类矩阵函数：（用来计算矩阵的 行列式、逆、秩、特征值 等等）
如,det inv rank eig poly 等等。
行列式 逆 秩 特征值 特征多项式例：B=,计算B的行列式、逆、秩、特征值、特征多项式。
（课堂上教师详细做。讲解poly：多项式系数（从高次到低次）构成的向量）
四．多项式例：用Matlab表示多项式
解,P=[1,-6,0,2,-8];
Px=poly2str(P,’x’)
执行结果为 Px=x^4-6*x^3+2*x-8
说明：（1）一个多项式的本质核心是它的各次项系数，把这些系数写成一个行向量（从高次到低次）称为该多项式的系数向量，系数向量完全决定了多项式，上面第一个命令就是把系数向量记作变量P；
（2）poly2str(P,’x’) 就是让机器根据系数向量P、以及符号’x’写出多项式。若把此句改为PPP=poly2str(P,’t’)，则执行结果为PPP=t^4-6*t^3+2*t-8 。
习题：用Matlab表示多项式
例：已知6次多项式G(x)的6个根为
-0.38 0.5 4 4.77 0.3+0.4i 0.3-0.4i
请写出G(x)的表达式。
解：R=[-0.38,0.5,4,4.77,0.3+0.4i,0.3-0.4i];
P=poly(R); （根据根R计算出系数向量）
P1=real(P); （因机器计算有误差，P带有虚部。根据数学知识，当虚根共轭成对出现时，必为实系数多项式。指令real(P)就是提取P的实部）
Gx=polt2str(P1,’x’)
习题：5次多项式F(x)的5个根为 -2 -1 0 1 2,请写出G(x)的表达式。
7次多项式H(x)的7个根为 3.2,-1.8,1,i,32i,请写出H(x)的表达式。
小结：指令poly(A)，若A为方阵，则产生A的特征多项式的系数；
若A为行向量，则产生以A为根的多项式的系数。
例：
解,命令,
再命令：c=conv(a,b); (产生系数向量)
Cx=polt2str(c,’x’)
完毕。
习题：=？
指令 roots(P) 的功能是：计算出以P为系数向量的多项式的根。
例：（求的根）
指令 polyval(P,A) 的功能是：计算出以P为系数向量的多项式在A的每个元素处的值。
例：
,求
解：命令为 P=[6,-8,0,3,-2,1];B=[2,-4,6;-3,5,-7];polyval(P,B) 完毕。
指令 polyvalm(P,A) 的功能是：A是方阵，按照矩阵运算规则计算出矩阵多项式的值。
例：
,求
解：命令为 P=[6,-8,0,3,-2,1];A=[2,-4;-3,5];polyvalm(P,A) 完毕。
五．Matlab工作区下的操作（初级）
清屏,clc
显示全部变量名,who（简显） 或 whos（详细显示）
显示某个变量AAA的内容,AAA 或 disp(AAA)
清除变量：（1）只清除变量AAA，则命令是 clear(AAA)
（2）清除机器中的全部变量，则命令是 clear
数据的输出格式,format 短型
format long 长型
format rat 分数型六．图形功能
1．二维图形 （命令词,plot ）
例1：已知x取值为 1,2,3,4,5,6
y取值为 0,0.5,0.7,1.1,0.9,0.2
请以x为横坐标、y为纵坐标，画出折线图。
解：以下三个方法中的任意一个都可以，其结果完全一样
（1）x=[1,2,3,4,5,6];y=[0,0.5,0.7,1.1,0.9,0.2];plot(x,y)
（2）x=1:6;y=[0,0.5,0.7,1.1,0.9,0.2];plot(x,y)
（3）y=[0,0.5,0.7,1.1,0.9,0.2];plot(y)
例2：画出正弦曲线y=sin(x)在区间上的图形。
练习： 画图一图多线：（在同一个坐标图中画两条以上的曲线）
法一：例1：x=0:pi/10:2*pi;y1=sin(x);y2=cos(x);plot(x,y1,x,y2)
例2：x=0:pi/10:2*pi;y1=sin(x);y2=cos(x);
y3=(sin(x)+cos(x))/2;plot(x,y1,x,y2,x,y3)
例3：在同一图中画四条线，横坐标范围是[0，2]，函数分别是
y= y= y= y
x1=0:pi/10:2*pi;y1=(1+sin(x))/2;x2=[0,2*pi];y2=x2/(2*pi);
y3=x1.*(2*pi-x1)/pi^2;y4=[0.6,0.6];
plot(x1,y1,x2,y2,x1,y3,x2,y4)
法二：重做例3：只需把plot(x1,y1,x2,y2,x1,y3,x2,y4)改成
hold on,plot(x1,y1),plot(x1,y3),plot(x2,y2),plot(x2,y4),hold off
习题：在同一图中画四条线，横坐标范围是[1，2+1]，函数分别是
,,,
加“坐标网格”：
法一：在 plot(……) 之后紧跟 grid
演示：plot(x1,y1,x2,y2,x1,y3,x2,y4),grid
或者hold on,plot(x1,y1),plot(x1,y3),plot(x2,y2),plot(x2,y4),grid,hold off
法二：事后补加（先画好没有网格的图，再补加网格）
演示：先命令plot(x1,y1),再命令grid
加“标记”：（1）给图加标题：title(‘******’)
（2）给坐标轴加标记：横轴：xlabel(‘******’)
纵轴：ylabel(‘******’)
（3）给曲线加标记：text(aaa,bbb,‘******’)
说明：这四个命令的用法与前面grid相同（紧跟plot(….)，或事后补加）；
红色符号是本质（包括：小括号、逗号、单引号），不可改变；
***** 是你准备加注的字符，由你自己决定；
你准备把“曲线标记”放在哪个位置？ 答：aaa是横坐标，bbb是纵坐标。
演示：略写，教师在课堂机器上演示。
对坐标系进行控制：
axis equal 使得横轴与纵轴的单位长度相同；
axis square 使得图呈现正方形；
axis(aax,bbx,aay,bby) 实现你的愿望：使得横坐标从aax到bbx、纵坐标从aay到bby,
演示：略写，教师在课堂机器上演示。
多幅图形：在同一个画面上建立多个坐标系，每个坐标系中各自画图、互不干扰关键词：subplot(m,n,p),把画面分割成个图形区域，p代表第几个图。
例如：借用前面的 x1,x2,y1,y2,y3,y4,
subplot(1,4,1),plot(x1,y1)
subplot(1,4,2),plot(x2,y2),grid
subplot(1,4,3),plot(x1,y3),grid
subplot(1,4,4),plot(x2,y4)
结果真难看！ 改
subplot(2,2,1),plot(x1,y1)
subplot(2,2,2),plot(x2,y2),grid
subplot(2,2,3),plot(x1,y3),grid
subplot(2,2,4),plot(x2,y4)
结果好看了！
习题：同一个画面上6幅图：
2．三维图形
先不动脑筋，看别人做的一个三维曲面图，感受一下Matlab作图效果。
例：画出曲面
解：x=-7.5:0.5:7.5;y=x; [X,Y]=meshgrid(x,y);
R=sqrt(X.^2+Y.^2)+eps;Z=sin(R)./R;
mesh(X,Y,Z)
（1）带网格的曲面图 (关键词：mesh )
例如：画曲面
x=-0.5:0.1:0.5;y=x;
[X,Y]=meshgrid(x,y); （组建三维图形的X,Y数组）
Z=sqrt(1-X.^2-Y.^2);
mesh(X,Y,Z) （画图）
（2）曲面的等高线图 (关键词：contour3 )
续上题：把 mesh(X,Y,Z) 改为 contour3(X,Y,Z,5)
(其中，5 是准备画5条等高线)
（3）空间三维曲线图 (关键词：plot3 )
例：螺旋线的参数方程为：
试在范围内，画出该曲线。
解：t=-20:pi/30:20;x=sin(t);y=cos(t);z=t;plot3(x,y,z)
习题：曲面，分别画网格图和等高线图。
曲线 自己适当取t的范围。
七．程序
1．编写、保存、打开、修改、执行程序编写新程序：击活菜单File，选New，单击M—file，进入编辑区，可编写程序。
例1：（演示：怎样编写、保存、打开、修改、执行程序？）
请按照规律1+2=3，2+3=5，3+5=8，5+8=13写出30个数作为一个数组a=[1,2,3,5,8,13,……],
简要说明：（1）保存新文件：File----Save,文件名必须字母打头；
（2）打开旧文件：File----Open----文件名，打开以后就可修改内容；
（3）执行：Debug----Run,
2．关系运算，逻辑运算
（1）关系运算：（比较两个数）
共6种运算符,< > <= >= == ~=
小于 大于 小于或等于 大于或等于 等于 不等于
练习,1<2 1>2 2<2 2<=2 1==2 1~=2
课堂判断,？？？？？？
Matlab规定：当关系成立时，结果是1；当关系不成立时，结果是0,
课堂做，且机器检验,(1<2)=? (1>2)=? (2<2)=?
(2<=2)=? (1==2)=? (1~=2)=?
习题,a=1:5 b=5:-1:1
c=a>=3,d=a==3,e=a<b,f=a==b
先手工算出 数组c,d,e,f,再上机器检验。
（2）逻辑运算,共有3种,& | ~
与(并且) 或 非(不是)
Matlab规定：(1) 当关系真时，结果是1；当关系假时，结果是0 ；
（2）若关系运算和逻辑运算混合了，则，关系运算优先，逻辑运算靠后，整体还叫关系表达式。
学习体会：这个符号“=”，它不叫“等号”，而叫“赋值号”。关系符号、逻辑符号都比它优先。
练习：a=1:8 b=8-a c=~(a>4) d=(a>=3)&(b<6)
则 a=
b=
c=
d=
3．选择结构（又称：条件语句）
格式一,if <关系表达式> 若 关系表达式 的值是1，则执行语句1；
<语句1> 若 关系表达式 的值是0，则不执行语句1，而去
end 执行end后面的语句。
格式二,if <关系表达式> 若 关系表达式的值是1，则执行语句1,再跳到end后面；
<语句1> 若 关系表达式的值是0，则不执行语句1，而执行语句2
else
<语句2>
end
格式三,if <关系表达式1> 若 关系表达式1的值是1，
<语句1> 则执行语句1,再跳到end后面；
elseif <关系表达式2> 若 式1假且式2真，
<语句2> 则执行语句2,再跳到end后面；
elseif <关系表达式3> 若 式1假、式2假且式3真，
<语句3> 则执行语句3,再跳到end后面；
…… …… …… 略。
end
4．循环结构格式一,for <循环参数> = <初值>,<步长>,<终值>
<语句1>
end
说明：对于每个参数值，都重复执行一次语句1,(课堂上，借用前面程序作解释)
例2：请输入、并显示一个20阶方阵，它的元素满足
解,笨法一：A=[1,2,3,4,…,19,20;2,3,4,…,21;……]
方法二：编程，程序清单如下：
for i=1:20
for j=1:20
A(i,j)=i+j-1;
end
end
A
( 循环结构中套循环！)
有时，程序中必需循环结构，但找不道合适的循环参数，或不知道循环次数，就用下面格式：
格式二,while <关系表达式1>
<语句1>
end
说明：每次执行完语句1，都要倒回去检验关系表达式1，只要关系表达式1的值是1，就再执行语句1。就这样重复执行直到关系表达式1的值是0为止。初学者慎用此结构，避免陷入死循环。
死循环例子：a=0;
while a>=0
a=a+1;
end
例3：记，求自然数m，使得
解,(课堂现编)
5．关于函数M----文件格式,function <被输出的变量名>=<文件名>（全部被引入的变量名）
注：变量名之间用逗号隔开例如，函数文件 jiecheng.m 的程序内容为：
function J=jiecheng (n)
J=1;
for k=1:n
J=J*k;
end
准备好该文件后，回到Matlab命令区。此时，你想计算几的阶乘？ 想算 3!，5!，12!，100!，只需分别下命令 jiecheng (3),jiecheng (5),jiecheng (12),jiecheng (100) 即可。
所以，文件jiecheng.m就是计算n阶乘的函数文件。
八．举例（用Matlab作计算）
小题1：随机产生一个满足“标准正态分布”的100维数组，记作a,再把数组a的这100个数的次序完全颠倒，得到一个新数组，记作b，输出b,
小题2：随机产生一个满足“(0,1)内均匀分布”的型矩阵，求出其最大元素及其所处的位置。
小题3：随机产生一个满足“(1,10)内均匀分布”的20维数组，用起泡法对20个数由小到大排序,即将相邻两个数比较,将小的调到前头.（不准用sort命令）。
小题4：编程求
小题5,一球从100米高度自由落下,每次落地后反跳回原高度的一半,再落下,求它在第10次落地时,共经过多少米?第10次反弹有多高?
大题1：
汽车头部的车灯，其反射面为一个旋转抛物面，方程为
。
灯丝（即：线光源）是一个直线段、长度为4、位于Y轴上且以原点为中点。 设点P(0,t,0)是线光源上任一点，则 设点M()是旋转抛物面任一点。从P到M的射线的反射线与平面 的交点记为N(25000,y,z)，计算公式为
点P在线光源上遍历，点M在抛物面上遍历，则，对应的点N在平面 上移动轨迹产生的区域就是反射光亮区。取t以步长0.2从-2到2，以步长1从-14.4到6.6。每固定一个，分别再遍历以步长从0到,编程画出反射光亮区。
大题2：
x
0
3
5
7
9
11
12
13
14
15
y
0
1.2
1.7
2.0
2.1
2.0
1.8
1.2
1.0
1.6
用Lagrange插值法，求当x分别取0,0.2,0.4,0.6,0.8,……,15时，对应的函数值y，并画出函数y=y(x)的图象.
提示：Lagrange插值公式
,
大题3,已知：，测试得下面数据
x,0 1153 2045 2800 3466 4068 4621 5135 5619 6152
y,0 20 40 60 80 100 120 140 160 183.5,
试用“最小二乘法”估计参数，并算出节点处的总误差.
大题4：（数值积分的龙贝格算法）计算，解法如下：
(k=1,2,3,……)
每当算出一个，就立即按下面公式计算
计算，输出格式为 format long,算到为止。顺便看看与准确值有多大误差？
大题5,求方程的最小正根。
要求：先观察确定有根区间[a,b]，区间宽度不大于0.2 ；
以下工作由程序来做：用“二分法”把有根区间宽度缩小到0.01，再启动“牛顿迭代公式”计算10次，输出近似根，输出格式为 format long
提示：0是整数？（笔误，改正：0是正数？）
方程f(x)=0,一个近似根记为，求根的牛顿迭代公式为
,
九．补充内容：
用Matlab软件求解微分方程
1．解析解
（1）一阶微分方程求的通解：dsolve('Dy=1+y^2','x')
求的通解：dsolve('Dy=1+x^2-y','x')
求的特解：dsolve('Dy=1+y^2',’y(0)=1’,'x')
（2）高阶微分方程求解其中，，命令为：
dsolve('x^2*D2y+x*Dy+(x^2-0.5^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x')
求的通解，命令为：
dsolve('D3y-2*Dy+y-4*x=0','x')
输出为：
ans=8+4*x+C1*exp(x)+C2*exp(-1/2*(5^(1/2)+1)*x)+C3*exp(1/2*(5^(1/2)-1)*x)
（3）一阶微分方程组求的通解：[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','x')
输出为,f =exp(3*x)*(cos(4*x)*C1+sin(4*x)*C2)
g =-exp(3*x)*(sin(4*x)*C1-cos(4*x)*C2)
若再加上初始条件，则求特解：
[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1','x')
输出为,f =exp(3*x)*sin(4*x)
g =exp(3*x)*cos(4*x)
2．数值解
（1）一阶微分方程
现以步长h=0.1用“4阶龙格—库塔公式”求数值解：
先建立“函数M—文件”：function f=eqs1(x,y)
f=y-2*x/y;
再命令,格式为：
[自变量,因变量]=ode45(‘函数文件名’,节点数组,初始值)
命令为,[x,y]=ode45('eqs1',0:0.1:1,1)
若还要画图，就继续命令,plot(x,y)
（2）一阶微分方程组
只须向量化，即可用前面方法：
function f=eqs2(x,y)
f=[cos(x)+2*y(1)-y(2);sin(x)-y(1)+2*y(2)];
将此函数文件，以文件名eqs2保存后，再下命令：
[x,y]=ode45('eqs2',0:0.1:1,[0.2;0.3])
（注：输出的y是矩阵，第i列为函数的数值解）
要画图，继续命令：hold on,plot(x,y(:,1)),plot(x,y(:,2)),hold off
（3）高阶微分方程先化成一阶微分方程组，再用前面方法。
上机练习：
准备：令，化成
用机器,函数文件eqs3内容？
命令？
画图？
用Matlab作最小二乘曲线拟合
1．用n次多项式作最小二乘拟合已知，要从（即：全体次数不高于n的多项式集合）中找一个，使得在节点处的总误差达到最小。
Matlab命令格式：系数数组=polyfit(节点数组，函数值数组，次数n)
例1：对函数C=C(t)测量得下面一组数据：
t,1 2 3 4 5 6 7 8 9
C：4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50
试分别用1次、2次、6次多项式作拟合，并画图显示拟合效果。
clear
hold on
x0=1:9;y0=[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50];
for i=1:9
plot(x0(i),y0(i),'+')
end
a1=polyfit(x0,y0,1),a2=polyfit(x0,y0,2),a6=polyfit(x0,y0,6)
x=0:0.1:10;
y1=polyval(a1,x);y2=polyval(a2,x);y6=polyval(a6,x);
plot(x,y1,x,y2,x,y6)
hold off
为了准确判断拟合效果，需计算“节点处的总误差”：（续前面程序）
wc1=sqrt(sum((polyval(a1,x0)-y0).^2))
wc2=sqrt(sum((polyval(a2,x0)-y0).^2))
wc6=sqrt(sum((polyval(a6,x0)-y0).^2))
2．用一般函数作最小二乘拟合已知，要用一个函数来近似代表，此函数中含有几个待定参数，现在的任务是：确定参数的值，使得在节点处的总误差达到最小。
用Matlab计算：先按照建立一个函数文件，文件第一句格式为：
function f=文件名(参数组,节点数组)
再用下面命令格式：
参数组=lsqcurvefit(‘函数文件名’，参数组初值，节点数组，函数值数组)
例2：对例1中的数据，用函数作拟合，并画图显示拟合效果，给出节点处总误差。
先建立并保存函数文件：文件名syp78hswj 内容为：
function f=syp78hswj(a,x0)
f=a(1)+a(2)*exp(-a(3)*x0);
再下面主程序：
clear
hold on
x0=1:9;y0=[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50];
for i=1:9
plot(x0(i),y0(i),'+')
end
cscz=[1,1,1];
a=lsqcurvefit('syp78hswj',cscz,x0,y0)
x=0:0.2:10;f=syp78hswj(a,x);plot(x,f)
f=syp78hswj(a,x0);
wc=sqrt(sum((f-y0).^2))
hold off
执行得：
节点处总误差 wc=0.0076
