matlab实验3-模型求解

题：用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
输入
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*= 不动点迭代：
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为矩阵，则显示每一列的正态 分布概率图形。
方差分析表中有6列： 第1列(source)显示：X中数据可变性的来源； 第2列(SS)显示：用于每一列的平方和； 第3列(df)显示：与每一种可变性来源有关的自由度； 第4列(MS)显示：是SS/df的比值； 第5列(F)显示：F统计量数值，它是MS的比率； 第6列显示：从F累积分布中得到的概率，当F增加时，p值减少。
A11x1 A12x2 A13x3 A21x1 A22x2 A23x3 An1x1 An2x2 An3x3
A1nxn B1 A2nxn B2
Annxn Bn
可以用矩阵形式表示为
即： AXB 则： XA\B
A11 A12 A13 A1n x1 B1 A21 A22 A23 A2n x2 B2
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近似值
多步法；Adams算法；高低精度均可到 10- 计算时间比 ode45 短 3～10-6
采用梯形算法
适度刚性情形
多步法；Gear’s 反向数值微分；精度中等 若 ode45 失效时，可尝试使用
单步法；2 阶Rosebrock 算法；低精度
当精度较低时，计算时间比 ode15s 短
梯形算法；低精度
当精度较低时，计算时间比 ode15s短
第三讲 ＭＡＴＬＡＢ模型求解 ——（常微分、偏微分）方程（组），概率统计
1、实验目的 2、实验内容 常微分方程求解函数 非线性方程求根 解（非）线性方程组 偏微分函数求解 概率统计函数 3、上机实践题 编程题 实践例题
内容提要
1、 实验目的
➢ 理解和运用MATLAB的常微分方程函数； ➢ 理解和运用MATLAB的非线性方程函数； ➢ 理解和运用MATLAB的（非）线性方程组函数； ➢ 理解和运用MATLAB的偏微分方程函数； ➢ 理解和运用MATLAB的概率统计函数； ➢ 理解和掌握简单处理程序设计
Baidu Nhomakorabea
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 ″,..]
函数 ode45
ode23
ode113
ode23t ode15s ode23s
ode23tb
ODE类型 非刚性
非刚性
非刚性
适度刚性 刚性 刚性
刚性
特点
说明
单步法；4，5 阶 R-K 方法；累计截断误差 大部分场合的首选方法 为 (△x)3
单步法；2，3 阶 R-K 方法；累计截断误差 使用于精度较低的情形 为 (△x)3
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)在每个迭代点的值
% 用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二阶导数函数图')
具体求解过程为：
Q 1Q 0f(Q 0,t0)h Q 2Q 1f(Q 1,t1)h Q3 Q2 f (Q2 , t2) h
简单欧拉方法程序
function [outx,outy]=MyEuler(fun,x0,xt,y0,PointNum) %MyEuler 用前向差分的欧拉方法解微分方程 %fun 表示f（x，y） %x0,xt表示自变量的初值和终值 %y0表示函数在x0处的值,其可以为向量形式 %PointNum表示自变量在[x0,xt]上取的点数 if nargin<5 | PointNum<=0 %如果函数仅输入4个参数值，则PointNum默认值为100
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页
h = normplot(X) %返回绘图直线的句柄
说明 样本数据在图中用“+”显示；如果数据来自正态分布，则图形显示为直线，而其它分布可 能在图中产生弯曲。
（3）绘制威布尔(Weibull)概率图形 函数 weibplot 格式 weibplot(X) %若X为向量，则显示威布尔(Weibull)概率图形，若X为矩阵，则显示每一列的威布尔 概率图形。 h = weibplot(X) %返回绘图直线的柄 说明 绘制威布尔(Weibull)概率图形的目的是用图解法估计来自威布尔分布的数据X，如果X是威布尔分布 数据，其图形是直线的，否则图形中可能产生弯曲。
f=f(:)'; y(k + 1,:) =y(k,:) +h*f; %对于所取的点x迭代计算y值
end
outy=y;
outx=x; %plot(x,y)%画出方程解的函数图
2、不动点迭代法求解非线性方程
迭代法是一种逐次逼近的方法，用某个固定公式反复校正根的近似值，使之逐步精确化，最后得 到满足精度要求的结果。
（4）样本的概率图形 函数 capaplot 格式 p = capaplot(data,specs) %data为所给样本数据，specs指定范围，p表示在指定范围内的概率。 说明 该函数返回来自于估计分布的随机变量落在指定范围内的概率
（5）附加有正态密度曲线的直方图 函数 histfit 格式 histfit(data) %data为向量，返回直方图 和正态曲线。 histfit(data,nbins) % nbins指定bar的个数， 缺省时为data中数据个数的平方根。
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函数产生两个图：标准的方差分析表图和盒图。
y1 y2
y3
0 0 0
1 0 3 x3
0 0
1 Y 0
2
x3
3e2x
x3
function dy=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的三阶导数
2.2 非线性方程求根函数
x = fzero(FUN,x0) %x0可以是数，或区间 x = fzero(FUN,x0,options) [x,fval]= fzero(FUN,x0,options) [x,fval,exitflag] = fzero(FUN,x0,options)
2.3 （非）线性方程组求解函数 线性代数方程组
2、 双因素方差分析 函数 anova2 格式 p = anova2(X,reps)
p = anova2(X,reps,'displayopt')
[p,table] = anova2(…)
[p,table,stats] = anova2(…) 说明 执行平衡的双因素试验的方差分析来比较X中两个或多个列（行）的均值，不同列的数据表示因素A的 差异，不同行的数据表示另一因素B的差异。如果行列对有多于一个的观察点，则变量reps指出每一单元观 察点的数目，每一单元包含reps行，如：
若对任意 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)