Matlab（9）
自动控制原理（3）
&[time,x]=solver(str,t,x0)&& 计算ODE或由字符串str 给定的ODE的值，部分解已在向量time中给出。在向量time
中给出部分解，包含的是时间值。还有部分解在矩阵x中给出，x的列向量每
个方程在这些值下的解。对于标量问题，方程的解将在向量x中给出。这些解
在时间区间t(1)到t(2)上计算得到。其初始值是x0 即x(t(1)).此方程组有str 指定的
M文件中函数表示出。这个函数需要两个参数：标量t 和向量x, 应该返回向量
x’(即x的导数)。因为对标量ODE来说，x和x’都是标量。在M文件中输入odefile
可得到更多信息。同时可以用命令numjac来计算Jacobi函数。
[t,x]=solver(str,t,x0,val)&& 此方程的求解过程同上，结构val包含用户给solver的命令。参见odeset和表1，可得到更多信息。
Ode45&&&&&此方法被推荐为首选方法。
Ode23&&&&&这是一个比ode45低阶的方法。
Ode113&&&&&用于更高阶或大的标量计算。
Ode23t&&&&&用于解决难度适中的问题。
Ode23s&&&&&用于解决难度较大的微分方程组。对于系统中存在常量矩阵的情况也有用。
Ode15s&&&&&与ode23相同，但要求的精度更高。
Ode23tb&&&&&用于解决难度较大的问题，对于系统中存在常量矩阵的情况也有用。
Set=odeset(set1,vak1,set2,val2,…)& 返回结构set,其中包含用于ODE求解方程的设置参数，&& 有关可用设置
的信息参见表1。
Odeget(set,’set1’)&&& 返回结构set中设置set1的值。
有许多设置对odeset控制的ODE解是有用的，参见表1。例如，如果要在求解过程中画出解的图形，可以
输入：inst=odeset(‘outputfcn’,’odeplot’);.也可使用命令odedemo。
表1&&& ODE求解方程的设置参数
RelTol&&&给出求解方程允许的相对误差
AbsTol&&给出求解方程允许的绝对误差
Refine&&给出与输入点数相乘的因子
OutputFcn&这是一个带有输入函数名的字符串，该字符串将在求解函数执行的每步被调用：odephas2(画出2D
的平面相位图)。Odephas3( 画出3D 的平面相位图)，odeplot( 画出解的图形)，odeprint( 显示中间结
OutputSel&是一个整型向量。指出哪些元素应该被传递给函数，特别是传递给OutputFcn
Stats&&如果参数Stats为on,则将统计并显示出计算过程中资源消耗情况
Jacobian… 如果编写ODE文件代码以便F（t,y,jocobian）返回dF/dy,则将jacobian设置为on
Jconstant…如果雅可比数df/dy是常量，则将此参数设置为on
Jpattern…& 如果编写ODE文件的编码以便函数F（[],[],jpattern）返回带有零的稀疏矩阵并输出非零元素dF/fy,则
需将Jpattern设置为on
Vectorized…如果编写ODE文件的编码以便函数F(t,[y1,y2……])返回[F(t,y1)& F(t,y2)…],则将此参数设置成on
Events…&& 如果ODE文件中带有参数‘events’，则将此参数设置成on
Mass…&& 如果编写ODE文件编码以实现函数F(t,[],‘mass’)返回M和M(t),应将此参数设置成on
MassConstant…如果矩阵M(t)是常量，则将此参数设置成on
MaxStep…&& 此参数是限定算法能使用的区间长度上限的标量
InitialStep…&& 给出初始步长的标量。如果给定的区间太大，算法就使用一个较小的步长
MaxOrder…&& 此参数只能被ode15s使用，它主要是指定ode15s的最高阶数,并且此参数应是从1到5的整数
BDF…&& 此参数只能被ode15s使用，如果倒推微分公式而不是用通常所使用的微分公式，则要将它设置为
NormControl…如果算法根据norm(e)&=max(Reltol*norm(y),Abstol)来步积分过程中的错误，则要将它设置成on
求解x’=-x^2 初值x(0)=1
创建函数xprim1,将此函数保存在M文件xprim1.m中：
function xprim=xprim1(t,x)
xprim=-x.^2;
然后调用MATLAB的ODE算法求解方程。然后画出解的图形：
[t,x]=ode45(‘xprim1’,[0,1],1);
plot(t,x,’-‘,t,x,’o’);
xlabel(‘time t0=0,tt=1’);
ylabel(‘x values x(0)=1’);
首先创建xprim2,将此函数保存在M文件xprim2.m中：
function xprim=xprim2(t,x)
xprim=x.^2;
然后调用MATLAB的ODE算法求解方程。然后画出解的图形：
[t,x]=ode45(‘xprim2’,[0,0.95],1);
注意：在MATLAB中计算出的点在微分绝对值大的区域内更密集些。
如果想在图像中直接观察，则在代码后输入以下命令
plot(t,x,’o‘,t,x,’-’);
xlabel(‘time t0=0,tt=0.95’);
ylabel(‘x values x(0)=1’);
例2求解下列二元一阶微分方程组
X1’=X1-0.1X1X2+0.01t
X2’=-X2+0.02X1X2+0.044
这个方程组用在人口动力学中。可以认为是单一化的捕食者---被捕食者模式。例如，狐狸和兔子。表示被捕食者， 表示捕食者。如果被捕食者有无限的食物，并且不会出现捕食者。于是有，这个式子是以指数形式增长的。大量的被捕食者将会使捕食者的数量增长；同样，越来越少的捕食者会使被捕食者的数量增长。而且，人口数量也会增长。
创建xprim3,将此函数保存在M文件xprim3.m 中：
function xprim=xprim3(t,x)
xprim=[x(1)-0.1*x(1)*x(2)+0.01*t; …
-x(2)+0.02*x(1)*x(2)+0.04*t];
然后调用一个ODE算法和画出解的图形：
[t,x]=ode45(‘xprim3’,[0,20],[30;20]);
plot(t,x);
xlabel(‘time t0=0,tt=20’);
ylabel(‘x values x1(0)=30,x2(0)=20’);
在MATLAB中,也可以根据x2函数绘制出x1图形,命令plot(x(:2),x(:1)) 可绘制出平面相位图.
例3 带参数的微分方程
X1’=a-(b+1)X1+X1^2X2
X2’=bX1+X1^2X2
方程由下面的M文件stiff1.m定义:
functionstiff=stiff1(t,x)
&& %变量不能放入参数表中
stiff=[0,0];&&& %stiff必须是一个冒号变量
stiff(1)=a-(b+1)*x1+x(1)^2*x(2);
stiff(2)=b*x(1)-x(1)^2*x(2);
下面的M文件给出一个比较困难的问题：
[t,X]=ode23(‘stiff1’，[0 10],[1 3]);
运行后得到的结果如下：
elapsed_time=72.1647
用专门解决复杂问题的解法ode23s,将得到较好的结果：
elapsed_time=1.0098
对于边界值问题，除了微分方程，还有边界处的值。在一维下这意味着至少有两个条件。
例4假设要研究一根杆的温度分布情况。这根杆一端的温度是T0,另一端的温度是T1。令y(x)表示这根杆的温度，函数f(x)表示加热源。从时间t=0开始，在相当长的时间内加热这根杆，直至达到平衡状态。这就是所谓的定常值或稳定状态。这个定常值可由下面的方程模型表示：
-Y’(X)=f(X),0&X&1
假设这根杆两端为：x=0和x=1。
假设在其两端又一根固定的柱子（或者可以看成是一个连接两个岛屿的桥），如图7所示。
令y(x)表示加载函数g(x)后弯曲的柱子。此问题需要有两个关于此柱子两端的边界条件。假设这根柱子非常牢
固的固定在墙上，即y在墙上的导数是0。可以得到下面的ODE,其中介绍了自然协调系统：
Y’’’(X)=g(x),0&X&1
Y(0)= Y’(0)= Y(1)= Y’(1)=0
解y(x)的导数可由有限的差分代替
如果用这些差分方程来代替ODE中的导数，就能得到一个所有未知的yi的方程组。其系数矩阵是一个有序区
间，此区间的宽度决定于这个微分方程的导数个数.
根据前面的温度模型的方程研究一下杆的温度分布，将所有的导数换成不同的差分并得到：
其中fj=f(xj)。为了简单起见，设M=6,即给定的y0和y6，而y1,y2,….y5 为未知变量。于是就有
注意，y0=T0和yM=T1必须移到方程组的右边。此时得到的矩阵是一个对角矩阵，其对角线上的元素为2，并且上一对角线和下一对角线上的元素为1。
下面解此问题的文件temperature.m。用户必须先给出分段数及f(x)(用点符号)，最后给出T0和T1。有关稀疏矩阵的更多信息参见其他资料。
%杆上的温度分布，用T0和T1分别表示两端温度
%这根杆放在x坐标的0和1 区间上，并被分成M个子区间，每个子区间的长度为1/M
%创建稀疏矩阵方程Ax=b并求解
%矩阵A对角阵，并以稀疏矩阵的形式存储
M=input(‘Give the number ofsubintervals (M):’);
Deltax=1/M;
xx=0:deltax:1;
funcStr= input(‘give f(x),the extra heatsource(e.g.,x.^3):’,’s’);
T0=input(‘Give y(0) (left): ’);
T1=input(‘Give y(1) (right): ’);
%构造对角阵和方程右边b
vectorOnes=ones(M-1,1);
A=spdiags([-vectorOnes,2*vectorOnes,-vectorOnes],[-1 0 1],M-1,M-1);
x=xx(2:end-1);
%x为区域内的值。
f=eval(funcStr);
%响应的f(x)的值。
b=deltax^2f;
b(1)=b(1)+T0; %对边界值x=0,x=1 进行特殊处理。
b(end)=b(end)+T1;
%解线性方程
%y在区间内：j=1:M-1.
y=[T0;y;T1];
%y在整个区间内：0&=x&=1.
%上面图形表示外部热源。
%下面图形表示杆上的热分布。
subplot(2,1,1);
plot(x,f);
title(‘External heatsource f(x).’,’FontSize’,14);
subplot(2,1,2);
plot(xx,y,’r’);
title(‘Tempearture distribution in a rod.’,’FontSize’,14);
%将区间分成等份，根据方程f(x)=在图中可以得到解。
例5 Matlab 中二阶常微分方程的数值解法
y''-3y'+2y=1;& y(0)=1;&y'(0)=0; 用数值解法求y(0.5)
odefun=@(t,x)[x(2);3*x(2)-2*x(1)+1];
[t,y]=ode45(odefun,[0:0.01:2],[1 0]);
y= dsolve('D2y-3*Dy+2*y=1','Dy(0)=0','y(0)=1');
y =exp(t) - exp(2*t)/2 + 1/2
&& feval(@(t)exp(t) - exp(2*t)/2 +1/2,0.5)
ans = 0.7896
例6求解范德普方程
首先是函数定义：
function xdot = myFcn (t, x)
xdot = [x(2) ; u(1-x(1)^2)*x(2)-x(1)]
其次是主程序：
% Solving options predetermined
options=odeset('OutputFcn',@odephas2,'reltol',1e-5);
% Solving vdp equation withdifferent ini conditions
for a = -3 : 0.1 : 3
b1=sqrt(25-a^2);
b2=-sqrt(25-a^2);
% Solve the problem
[t,y]=ode45(@myFcn, [0 20],[a b1], options);
[t,y]=ode45(@myFcn, [0 20],[a b2], options);
例7x1.^2-10x1+x2.^2+b=0
&&&& x1*x2.^2+x1-10x2+8=0
用MATLAB实现，编写非线性方程组的M文件fx1.m如下所示：
&& function y=fx1(x)
&&y(1)=x(1)*x(1)-10*x(1)+x(2)*x(2)+8;
&&y(2)=x(1)*x(2)*x(2)+x(1)-10*x(2)+8;
编写非线性方程组导数的M文件dfx1.m
& function y=dfx1(x)
&&&& y(1)=2*x(1)-10;
&&&& y(2)=2*x(2);
&&&& y(3)=x(2)*x(2)+1;
&&&& y(4)=2*x(1)*x(2)-10;
附一篇文章《用MATLAB求解非线性微分方程》
原始地址：http://blog.csdn.net/anhuixue/article/details/7560558
总结一下MATLAB中求解微分方程的思路和步骤。MATLAB中的帮助文件，内容实在是特别给力，简明的不能再简明了。可惜了是英文写的，即使是天天要面对英语的我，也颇为头疼。那么以下我将结合我对MATLAB中帮助文档的理解，提出我对求解最简单微分方程的方法的总结。以求解一个经典的非线性方程为例。
一求解一下这个方程，判断她是否稳定。要是稳定，那么她是否存在极限环;
二.一般的计算机求解方程的方法：
首先把该方程改写成一个规范的形式，一般使用状态空间表示法；
调用已有的算法进行求解；
最后对得出的结果进行处理，比如画图之类的。
三.输入待求解的方程。
MATLAB中关于输入输入待解方程的语句特别简单。需要先定义一个普通函数，函数名任意，姑且叫做myFcn，格式如下
function xdot = myFcn (t, x)&
需要注意的是，函数必须含有t, x两个参数，名称可以自己任意定。xdot是这个函数本身的返回值，只出现在这个函数内部，因此也可以任意定。
定义中的x被视为是一个列向量，x(i)表示列向量中的第i个分量。那么F函数的每一个分量即简单地用表达时给出即可。其中的自变量可以引用x(i)。以范德普方程为例：
xdot = [x(2) ; u(1-x(1)^2)*x(2)-x(1)]&
于是，这两句话便构成了待解函数。
四.调用MATLAB函数进行求解
通常人工求解微分方程需要知道初始值，计算机求解也不例外。另外，由于非线性方程一般只有数值解，故计算精度也可以调整。这些都是可以自己调整的参数。
调用MATLAB计算求解常微分方程的模式很简单，格式为：
[t, x] = ode45(@myFcn, [开始时间 结束时间], [初始值列向量],options)
注意到求解出来的结果是[t, x]即一堆数，所以需要我们进行后处理比如画图之类的。
ode45表示了MATLAB中的一种内置计算微分方程的算法，最为常用，出于娱乐目的，就先用这个。
@表示待求解的方程，如myFcn。
[开始时间 结束时间]表示求解的时间段，不必定义间隔。
[初始值列向量]就不用多说了。通常在求解之前定义一个变量x0列向量表示初始值，然后输入起来就方便多了。
options本身是一个变量。注意，她的名称不固定，但是他是一个以结构体为类型的变量。options很重要，稍后讨论。
五.进行后处理。在得到[t, x]后可以自己用plot神马的进行画图。
六.options的用法。
options在MATLAB帮助中是这样定义格式的：
&&&&&&options = odeset (“Name”, Value, “Name”, Value, …)
意思是，options这个变量要用到odeset这个函数来对他进行赋值。而odeset这个函数的参数必须是成对出现的：一个名称对应一个数值。
我们要用到的是这样的：
options= odeset(“OutputFcn”, @odephas2,“reltol”, 1e-5);
注意，odeset中的参数名称不可任意定，因为都是预定义好的。”OutputFcn”使用预定义的函数进行输出，而与定义的函数有好几种，使用@进行选择，我们选odephas2即画相平面图(phase portrait)。之后的那个参数表示相对误差的最大允许值，不说也明白。
有一个问题，用odephas2画出的图图不好看，因为上面打了好多圈圈。解决办法是找到该文件，修改其中所有的plot语句即可，把那个”-o”去掉就行了。
参考知识库
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：14333次
排名：千里之外
原创：24篇
转载：13篇
(1)(6)(9)(14)(7)数值分析_第三次上机_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
数值分析_第三次上机
上传于||暂无简介
阅读已结束，如果下载本文需要使用0下载券
想免费下载更多文档？
定制HR最喜欢的简历
下载文档到电脑，查找使用更方便
还剩2页未读，继续阅读
定制HR最喜欢的简历
你可能喜欢Custom x-axis values in a matlab plot - Stack Overflow
to customize your list.
Announcing Stack Overflow Documentation
We started with Q&A. Technical documentation is next, and we need your help.
Whether you're a beginner or an experienced developer, you can contribute.
Currently when I plot a 9 by 6 array, the x-axis of the figure is just 1, 2, 3 up to 9. The Y-axis shows the correct values.
Instead of 1 to 9 I would like the x-axis values to be custom. They should be
100 200 400 00
instead. I tried
set(gca,'XTick', [100 200 400 00
But that's not the correct way to do it. Is there a Matlab option to have these custom values for the x-axis? Why is Matlab just using 1 to 9 anyway?
7,33011223
If you want to keep distances between x-values (e.g. 1:9) and only change the labels (not the distances between x-values), try this:
y = rand(9,6);
labels = [100 200 400 00
set(gca, 'XTick', 1:length(labels)); % Change x-axis ticks
set(gca, 'XTickLabel', labels); % Change x-axis ticks labels.
You should be using
instead of XTick.
MATLAB plots every column as a seperate curve. So, that means you have 6 curves and 9 data points for each curve. x-axis data is 1-9 because you did not provide any data for MATLAB to plot with.
Furthermore, you probably want the wrong thing. Doing this will give you equal spacing. It will just replace 1-9 with your array. Since your x-axis data is not equally spaced, it will be weird.
You may want to do it like this:
xdat = [100 200 400 00
ydat = rand(9,6); % Your y-axis data
plot(xdat, ydat)
x = [100 200 400 00
y = rand(9,6); % Your y-axis data
plot(x, y);
set(gca,'XTick',x); % Change x-axis ticks
set(gca,'XTickLabel',x); % Change x-axis ticks labels to desired values.
Please note that due to the very different values/magnitudes to use un x-axis you can get some x-labels very close (and unreadable)
Your Answer
Sign up or
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Post as a guest
By posting your answer, you agree to the
Not the answer you're looking for?
Browse other questions tagged
Stack Overflow works best with JavaScript enabled利用MATLAB求解这样一个非线性自回归模型，y=ax^2+bx+c ,求系数a,b,c的值，感激不尽_百度知道
利用MATLAB求解这样一个非线性自回归模型，y=ax^2+bx+c ,求系数a,b,c的值，感激不尽
求a其中x=[100 105 96 102]
x^2=[ ]y=[105 96 102 101]不知道该怎么做,b，谢谢大家了。开始我用曲线拟合的方法做出来结果不对，不知道该怎么做了,c
不用二次拟合的方法能做吗？
提问者采纳
.5665c = -2： a =
-0;X=[ones(size(x))t=regress(y.2353b =
46.^2+b*x+c----------------------------------结果如下.^2]; x=[100 105 96 102]'a=t(3)b=t(2)c=t(1) a*x.7066
102;y=[105 96 102 101]&#39.0099
95;&gt.1001&&gt.ans =
提问者评价
其他类似问题
为您推荐：
其他2条回答
a(1)*x^2+a(2)*x+a(3)'a'x'a=lsqcurvefit(,x;& a=a(1)a =
-0; y=[105 96 102 101];;&gt,b=a(2),');c=a(3), 1],y);fun=inline(&#39.5663&gt, 1,&#39x=[100 105 96 102] ,a=a(1)c = -2; b=a(2)b =
多项式拟合 x=[100 105 96 102];y=[105 96 102 101];t=polyfit(x,y,2);a=t(1)b=t(2)c=(3)
自回归模型的相关知识
等待您来回答
下载知道APP
随时随地咨询
出门在外也不愁用matlab求解常微分方程_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
用matlab求解常微分方程
上传于||文档简介
&&基​础​性​的​m​a​t​l​a​b​了​解​资​料​,​初​步​熟​悉​m​a​t​l​a​b​解​微​分​方​程​的​方​法
阅读已结束，如果下载本文需要使用0下载券
想免费下载更多文档？
定制HR最喜欢的简历
下载文档到电脑，查找使用更方便
还剩2页未读，继续阅读
定制HR最喜欢的简历
你可能喜欢