matlab模型求解

2.3 （非）线性方程组求解函数 线性代数方程组
A11 x1 A12 x2 A13 x3 A1n xn B1 A x A x A x A x B 21 1 22 2 23 3 2n n 2 An1 x1 An 2 x2 An 3 x3 Ann xn Bn
格式 p = capaplot(data,specs) %data为所给样本数据，specs指定范围， p表示在指定范围内的概率。
说明 该函数返回来自于估计分布的随机变量落在指定范围内的概率
（5）附加有正态密度曲线的直方图
函数 histfit
格式 histfit(data) 和正态曲线。 %data为向量，返回直方图
格式 [x,y]=ode45(′ fun′, [x0,xn], y0,option]
说明： ① 适用于求解一阶常微分方程组 ② fun定义微分方程组的函数文件名 ③ [x0,xn]求解区域 ④ y0初始条件向量 ⑤ option可选参数，由ODESET函数设置，比较复杂 ⑥ x输出自变量向量，y输出[y, y ′, y ″,..]

1、 实验目的
理解和运用MATLAB的常微分方程函数；
理解和运用MATLAB的非线性方程函数； 理解和运用MATLAB的（非）线性方程组函数； 理解和运用MATLAB的偏微分方程函数；
理解和运用MATLAB的概率统计函数；
理解和掌握简单处理程序设计
2、实验内容
2.1 常微分方程函数 ode45 ode23 ode113 ode15s ode23s ode23t ode23tb
输出
迭代N 0 次还没有达到
精度要求信息
ຫໍສະໝຸດ Baidu
将连续函数方程f(x)＝0改写为等价形式：x=(x) 其中(x)也是连续函数，称为迭代函数。 不动点：若x*满足f(x*)=0，则x*=(x*)；反之，若 x*=(x*) ，则f(x*)=0 ，称x*为(x)的一个不动点。 不动点迭代： xk 1 ( xk ) (k=0,1,……) 若对任意 x0[a,b]，由上述迭代得序列{xk}，有极限
适度刚性 采用梯形算法 刚性
ode23s
ode23tb
刚性
刚性
单步法；2 阶Rosebrock 算法； 当精度较低时，计算时 低精度 间比 ode15s 短
梯形算法；低精度 当精度较低时，计算时 间比ode15s短
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、不动点迭代法求解非线性方程 迭代法是一种逐次逼近的方法，用某个固定公 式反复校正根的近似值，使之逐步精确化，最 后得到满足精度要求的结果。
输入
x1 , , N 0
k 1,2,, N 0
计算
xk 1 x g x1
是 输出
x x1
否
k, x
x1 x k N0
函数
ODE类 型
非刚性 非刚性 非刚性
特点
单步法；4，5 阶 R-K 方法； 累计截断误差为 (△x)3 单步法；2，3 阶 R-K 方法； 累计截断误差为 (△x)3
说明
大部分场合的首选方法 使用于精度较低的情形
ode45 ode23
ode113 ode23t ode15s
多步法；Adams算法；高低精 计算时间比 ode45 短 度均可到 10-3～10-6 适度刚性情形 多步法；Gear’s 反向数值微分； 若 ode45 失效时，可尝 精度中等 试使用
y1 0 dY d y2 0 dx dx y3 0
1 0 3 x3
0 0 1 Y 0 3e2 x 2 3 3 x x
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的三阶导数 % 用ode23 ode45 ode113解多阶微分方程 clear,clc
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函数产生两个图：标准的方差分析表图和 盒图。
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页
题：用MATLAB函数ode23,ode45,ode113,求解多阶常微分方程：
作如下近似：
Qn Q(tn )
则得到欧拉解法递推公式的一般形式：
Qn1 Qn f (Qn , tn ) h
具体求解过程为：
Q1 Q0 f (Q0 , t0 ) h Q2 Q1 f (Q1 , t1 ) h
Q3 Q2 f (Q2 , t2 ) h
3 2 d y d y dy 3 2x x 2 3 3 e 3 2 dx dx dx y (1) 1, y '(1) 10, y "(1) 30, x [1,10]
dy1 dx y2 dy2 y3 dx dy3 2 3 3e 2 x 3 y3 3 y2 3 x x dx x
方差分析表中有6列： 第1列(source)显示：X中数据可变性的来源； 第2列(SS)显示：用于每一列的平方和； 第3列(df)显示：与每一种可变性来源有关的自由度； 第4列(MS)显示：是SS/df的比值； 第5列(F)显示：F统计量数值，它是MS的比率； 第6列显示：从F累积分布中得到的概率，当F增加时， p值减少。
2、 双因素方差分析 函数 anova2 格式 p = anova2(X,reps) p = anova2(X,reps,'displayopt')
[p,table] = anova2(…)
[p,table,stats] = anova2(…) 说明 执行平衡的双因素试验的方差分析来比较X中两个或多个 列（行）的均值，不同列的数据表示因素A的差异，不同行的 数据表示另一因素B的差异。如果行列对有多于一个的观察点， 则变量reps指出每一单元观察点的数目，每一单元包含reps行， 如： reps=2 其余参数与单因素方差分析参数相似。
lim xk x *
k
则称迭代过程收敛，且x*=(x*)为(x)的不动点。
function [root,n]=StablePoint(f,x0,eps) if(nargin==2) 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、概率统计图
（1） 最小二乘拟合直线 函数 lsline 格式 lsline %最小二乘拟合直线 h = lsline %h为直线的句柄
（2） 绘制正态分布概率图形
函数 normplot
格式 normplot(X) %若X为向量，则显示正态分布概率图形，若X 为矩阵，则显示每一列的正态分布概率图形。 h = normplot(X) %返回绘图直线的句柄 说明 样本数据在图中用“+”显示；如果数据来自正态分布，则图 形显示为直线，而其它分布可能在图中产生弯曲。
简单欧拉方法程序
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 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)在每个迭代点的值 f=f(:)'; y(k + 1,:) =y(k,:) +h*f; %对于所取的点x迭代计算y值 end outy=y; outx=x; %plot(x,y)%画出方程解的函数图
对于等间隔节点
t tn1 tn h tn1 tn h
可以得到： tn Q精确值 t0 t1 t2
n=0,1,2
…. …. ….
tn
…. …. ….
Q(t0) Q(t1) Q(t2)
Q(tn)
Q近似值
Q0
Q1
Q2
Qn
在tn节点上，微分方程可以写为
Q(tn1 ) Q(tn ) f Q(tn ) , tn h
（3）绘制威布尔(Weibull)概率图形 函数 weibplot 格式 weibplot(X) %若X为向量，则显示威布尔(Weibull)概率图形，若X 为矩阵，则显示每一列的威布尔概率图形。 h = weibplot(X) %返回绘图直线的柄 说明 绘制威布尔(Weibull)概率图形的目的是用图解法估计来自威布尔分 布的数据X，如果X是威布尔分布数据，其图形是直线的，否则图形中可 能产生弯曲。 （4）样本的概率图形 函数 capaplot
科学计算与MATLAB 实验讲义
中南大学材料科学与工程学院
第三讲 ＭＡＴＬＡＢ模型求解
——（常微分、偏微分）方程（组），概率统计
内容提要
1、实验目的 2、实验内容 常微分方程求解函数 非线性方程求根 解（非）线性方程组 偏微分函数求解 概率统计函数 3、上机实践题 编程题 实践例题
可以用矩阵形式表示为 即： 则：
格式
A11 A12 A21 A22 An1 An 2
A13 A23 An 3
A1n x1 A2 n x2 Ann xn
B1 B2 Bn
AX B
X A\ B
solve('eqn1','eqn2',...,'eqnN','var1,var2,...,varN')
histfit(data,nbins) % nbins指定bar的个数，
缺省时为data中数据个数的平方根。
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