f=f(:)'; y(k + 1,:) =y(k,:) +h*f; %对于所取的点x迭代计算y值
end
outy=y;
outx=x; %plot(x,y)%画出方程解的函数图
2、不动点迭代法求解非线性方程
迭代法是一种逐次逼近的方法，用某个固定公式反复校正根的近似值，使之逐步精确化，最后得 到满足精度要求的结果。
2、实验内容
2.1 常微分方程函数
ode45 ode23 ode113 ode15s ode23s ode23t ode23tb
格式 [x,y]=ode45(′ fun′, [x0,xn], y0,option]
说明： 适用于求解一阶常微分方程组 fun定义微分方程组的函数文件名 [x0,xn]求解区域 y0初始条件向量 option可选参数，由ODESET函数设置，比较复杂 x输出自变量向量，y输出[y, y ′, y ″,..]
% 用ode23 ode45 ode113解多阶微分方程 clear,clc [x23,y23]=ode23('myfun03',[1,10],[1 10 30]); [x45,y45]=ode45('myfun03',[1,10],[1 10 30]); [x113,y113]=ode113('myfun03',[1,10],[1 10 30]); 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二阶导数函数图')
若对任意 x0 [a,b]，由上述迭代得序列{xk}，有极限
(x*) ，则f(x*)=0 ，称x*为 (k=0,1,……)
(x)的一个不动点。
则称迭代过程收敛，且x*= (x*)为 (x)的不动点。
xk1(xk)
lk imxk x*
function [root,n]=StablePoint(f,x0,eps) if(nargin==2)
PointNum=100;
end if nargin<4 %y0默认值为0
y0=0;
end h=(xt-x0)/PointNum;%计算步长h x=x0+[0:PointNum]'*h;%自变量数组 y(1,:) = y0(:)';%将输入存为行向量，输入为列向量形式
for k = 1:PointNum f=feval(fun,x(k),y(k,:));%计算f(x,y)在每个迭代点的值
3、 上机实践
3.1编程题 1、欧拉数值算法（差分法）求解常微分方程 差分法就是用差商近似代替微商，即
Q dQ t dt
代入微分方程得到：
Q(t t)Q(t) f (Q,t) t
Q(t t) Q(t) f (Q,t)t
对于等间隔节点
ttn1tnh tn1 tnh
可以得到：
n=0,1,2
tn Q精确值 Q近似值
题：用MATLAB函数ode23,ode45,ode113,求解多阶常微分方程：
x3 d3y2d2y3dy3e2x dx3 dx2 dx
y(1)1,y'(1)10,y"(1)30,x[1,10]
dy1
dx
y2
dy2 dx
y3
dy3 dx
2 x3
y3
3 x3
y2
3e2x x3
dY dx
d dx
多步法；Adams算法；高低精度均可到 10- 计算时间比 ode45 短法；Gear’s 反向数值微分；精度中等 若 ode45 失效时，可尝试使用
单步法；2 阶Rosebrock 算法；低精度
当精度较低时，计算时间比 ode15s 短
梯形算法；低精度
当精度较低时，计算时间比 ode15s短
t0
t1
t2 …. tn ….
Q(t0) Q(t1) Q(t2) …. Q(tn) ….
Q0
Q1 Q2 …. Qn ….
在tn节点上，微分方程可以写为
Q ( tn 1 ) Q ( tn ) fQ ( tn ),tn h
作如下近似：
Qn Q(tn)
则得到欧拉解法递推公式的一般形式：
Q n 1Q nf(Q n,tn)h
reps=2 其余参数与单因素方差分析参数相似。
3、概率统计图
（1） 最小二乘拟合直线 函数 lsline 格式 lsline %最小二乘拟合直线
h = lsline %h为直线的句柄
（2） 绘制正态分布概率图形
函数 normplot
格式 normplot(X) %若X为向量，则显示正态分布概率图形，若X为矩阵，则显示每一列的正态 分布概率图形。
h = normplot(X) %返回绘图直线的句柄
说明 样本数据在图中用“+”显示；如果数据来自正态分布，则图形显示为直线，而其它分布可 能在图中产生弯曲。
（3）绘制威布尔(Weibull)概率图形 函数 weibplot 格式 weibplot(X) %若X为向量，则显示威布尔(Weibull)概率图形，若X为矩阵，则显示每一列的威布尔 概率图形。 h = weibplot(X) %返回绘图直线的柄 说明 绘制威布尔(Weibull)概率图形的目的是用图解法估计来自威布尔分布的数据X，如果X是威布尔分布 数据，其图形是直线的，否则图形中可能产生弯曲。
2、 双因素方差分析 函数 anova2 格式 p = anova2(X,reps)
p = anova2(X,reps,'displayopt')
[p,table] = anova2(…)
[p,table,stats] = anova2(…) 说明 执行平衡的双因素试验的方差分析来比较X中两个或多个列（行）的均值，不同列的数据表示因素A的 差异，不同行的数据表示另一因素B的差异。如果行列对有多于一个的观察点，则变量reps指出每一单元观 察点的数目，每一单元包含reps行，如：
函数 ode45
ode23
ode113
ode23t ode15s ode23s
ode23tb
ODE类型 非刚性
非刚性
非刚性
适度刚性 刚性 刚性
刚性
特点
说明
单步法；4，5 阶 R-K 方法；累计截断误差 大部分场合的首选方法 为 (△x)3
单步法；2，3 阶 R-K 方法；累计截断误差 使用于精度较低的情形 为 (△x)3
An1 An2 An3 Ann xn Bn
格式 solve('eqn1','eqn2',...,'eqnN','var1,var2,...,varN')
Matlab非线性方程组求解
格式 X=fsolve(FUN,X0)
2.4 概率统计函数
单因素方差分析 函数 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函数产生两个图：标准的方差分析表图和盒图。
输入
x1,, N0
k1,2,,N0
计算 xk1 x gx1
kN0
xx1
否
x1x
是 输出 k , x
输出
迭代 N0次还没有达到
精度要求信息
将连续函数方程f(x)＝0改写为等价形式：x= (x) 其中 (x)也是连续函数，称为迭代函数。
不动点：若x*满足f(x*)=0，则x*= (x*)；反之，若x*= 不动点迭代：
（4）样本的概率图形 函数 capaplot 格式 p = capaplot(data,specs) %data为所给样本数据，specs指定范围，p表示在指定范围内的概率。 说明 该函数返回来自于估计分布的随机变量落在指定范围内的概率
（5）附加有正态密度曲线的直方图 函数 histfit 格式 histfit(data) %data为向量，返回直方图 和正态曲线。 histfit(data,nbins) % nbins指定bar的个数， 缺省时为data中数据个数的平方根。
eps=1.0e-4; end tol=1; root=x0; n=0; while(tol>eps)
n=n+1; r1=root; root=subs(sym(f),findsym(sym(f)),r1)+r1; tol=abs(root-r1); end
3.2上机例题
1、常微分方程函数 275页，例题7-47 273-274页，例题7-45/46 2、（非）线性方程组求解例题 246页，例7-17/18 265页，例7-39 3、偏微分方程求解函数 338页，例题 4、概率统计函数 307页，例9-21；308页，例9-22；309页，例9-23 310页