科学计算：Python&VS.&MATLAB（5）----常微分方程数值解
科学计算：Python VS.
MATLAB（5）----常微分方程数值解
一、常微分方程的一般理论
凡含有参数，未知函数和未知函数导数 (或微分)
的方程，称为微分方程，有时简称为方程，未知函数是一元函数的微分方程称作常微分方程，未知数是多元函数的微分方程称作偏微分方程。微分方程中出现的未知函数最高阶导数的阶数，称为微分方程的阶。
二、使用Python
scipy中提供了用于解常微分方程的函数odeint()，完整的调用形式如下：
scipy.integrate.odeint(func,
y0, t, args=(), Dfun=None, col_deriv=0,
full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None,
h0=0.0, hmax=0.0,hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12,
mxords=5, printmessg=0)
实际使用中，还是主要使用前三个参数，即微分方程的描写函数、初值和需要求解函数值对应的的时间点。接收数组形式。这个函数，要求微分方程必须化为标准形式，即dy/dt=f(y,t,)。而我们会发现，高阶方程的标准化工作，其实是解微分方程最主要的工作。下面是一个二阶常微分方程的例子，通过这个例子，学会变化标准形式的方法。
典型的二阶方程：y"+a*y'+b*y=0
，不满足标准式。变量替换，假设输入是Y=[y,y’]，输出是Y’=[y’,y”]，那么方程就化为了Y’=F(Y,t)的形式。我们知道Y[0]=y,
Y[1]=y’，所以输出形式就是：
y”=-a*y’-b*y=-a*Y[1]-B*Y[0]
有了这个形式，我们很容易写出Python代码：
from scipy.integrate import odeint
from pylab import *
deriv(y,t):&&&&&&&
# 返回值是y和y的导数组成的数组
&&& return
array([ y[1], a*y[0]+b*y[1] ])
time = linspace(0.0,50.0,1000)
array([0.])&&&&
y = odeint(deriv,yinit,time)
plot(time,y[:,0],label='y')&&&
#y[:,0]即返回值的第一列，是y的值。label是为了显示legend用的。
hold('on')
plot(time,y[:,1],label="y'")&&&&
#y[:,1]即返回值的第二列，是y’的值
xlabel('t')
ylabel('y')
-------------------------------------------------------
三、使用MATLAB
同样，MATLAB解微分方程的第一步工作也是化简微分方程为一阶的标准形式，这个和Python是完全一样的。只是函数的编写不一样。MATLAB可以在m文件函数中定义一个函数和多个子函数，但是子函数只能由同一m文件中的函数调用。并且还要求m文件名必须和主函数同名；而Python无此限制，在此可见一斑，也能看出Python代码的相对更优美。另外，Python的索引从0开始，MATLAB的从1开始，为了这个习惯，编写函数的时候也要有所注意。
MATLAB在这里也有一个优势，就是提供了多个求解函数，以应对时间需求、精度需求和刚性问题等要求。这就比scipy中的函数灵活了不少。如果使用scipy，只能自己编写函数了！
这些函数几乎统一的调用方式是：[to,yo]=solver(func,tspan,y0,options)
其中，solver是ode45(最常用)、ode23、ode23t、ode23s、ode23tb等
例1：y”+4*y=3*sin(t)&
y(0)=1,y’(0)=0
定义一个函数，描写这个微分方程：
function dy=fun1(t,y)
dy(1,1)=y(2)
dy(2,1)=3*sin(t)-4*y(1)
将这个函数存为fun1.m,
然后写一个脚本调用它，解方程。不妨这个脚本就叫做callfun1.m
tspan=[0,5];
[t,y]=ode45(@fun1,tspan,y0);
plot(t,y(:,1),'k-');
plot(t,y(:,2),'k:');
xlabel('t');
画出来的图是这样子的：
下面用Python实现MATLAB部分的例子。
from scipy.integrate import odeint
from pylab import *
def deriv(y,t):
&&& return
array([ y[1], 3*sin(t)-4*y[0] ])
#设置初值和t区间
time = linspace(0.0,5.0,1000)
yinit = array([1,0])
#解方程并绘图
y = odeint(deriv,yinit,time)
plot(time,y[:,0],label='y')&
hold('on')
plot(time,y[:,1],label="y'")&
xlabel('t')
ylabel('y')
结果如图，可以和MATLAB绘制出来的对比一下：
-------------------------
下面用MATLAB实现Python部分的例子。
函数脚本（fun1.m）：
function dy=fun1(t,y)
dy(1,1)=y(2)
dy(2,1)=-2.0*y(1)-0.1*y(2)
调用脚本（callfun1.m）：
tspan=[0,50];
[t,y]=ode45(@fun1,tspan,y0);
plot(t,y(:,1),'k-');
plot(t,y(:,2),'k:');
xlabel('t');
绘图结果：
补充说明：
MATLAB还可以使用inline的形式调用函数。比如：
tspan=[0,5];
fun1=inline('[y(2);3*sin(t)-4*y(1)]','t','y');&&&
%注意这一句
[t,y]=ode45(fun1,tspan,y0);
plot(t,y(:,1),'k-');
plot(t,y(:,2),'k:');
xlabel('t');
这时候，不需要关心文件名的问题了，但是这种写法明显不利于推广和模块化，所以不推荐使用这样的形式。
我的更多文章：
( 19:46:30)( 13:03:50)( 13:41:19)( 16:45:22)
已投稿到：
以上网友发言只代表其个人观点，不代表新浪网的观点或立场。【图文】常微分方程求解_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
常微分方程求解
大小：1.27MB
登录百度文库，专享文档复制特权，财富值每天免费拿！
你可能喜欢温馨提示！由于新浪微博认证机制调整，您的新浪微博帐号绑定已过期，请重新绑定！&&|&&
LOFTER精选
网易考拉推荐
用微信&&“扫一扫”
将文章分享到朋友圈。
用易信&&“扫一扫”
将文章分享到朋友圈。
函数 、ode23、ode113、ode15s、ode23s、ode23t、ode23tb
功能 微分方程（）组初值问题的数值解
%在区间上，从到，用初始条件求解显式微分方程。对于标量与列向量，函数必须返回一的列向量。解矩阵中的每一行对应于返回的时间列向量中的一个时间点。要获得问题在其他指定时间点上的解，则令。
%用参数（用命令生成）设置的属性（代替了缺省的积分参数），再进行操作。常用的属性包括相对误差值（缺省值为）与绝对误差向量（缺省值为）。
& %将参数等传递给函数，再进行计算。若没有参数设置，则令。
参数说明：
为命令、。
为显式常微分方程，或为包含一混合矩阵的方程。命令只能求解常数混合矩阵的问题；命令与可以求解奇异矩阵的问题。
积分区间（即求解区间）的向量。要获得问题在其他指定时间点上的解，则令。
包含初始条件的向量。
用命令设置的可选积分参数。
传递给函数的可选参数。
．求解具体的具体步骤：
& （）根据问题所属学科中的规律、定律、公式，用微分方程与初始条件进行描述。
&&（）运用数学中的变量替换：，把高阶（大于阶）的方程（组）写成一阶微分方程组。
& （）根据（）的结果，编写能计算导数的函数文件。
& （）将m文件与初始条件传递给求解器中的一个，运行后就可得到的、在指定时间区间上的解列向量
．求解器见下列表示（表来源于网络）：
& 下面我们做应用举例：
& clear all close all clc t_final=100; x0=[0;0;1e-10]; [t,x]=ode45('lorenzeq',[0,t_final],x0); plot(t,x),
plot(t,x(:,1)),
plot(t,x(:,2)),
plot(t,x(:,3)),
plot3(x(:,1),x(:,2),x(:,3)); axis([10 40 -20 20 -20 20]);
& %%lorenzeq.m function xdot=lorenzeq(t,x) xdot=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)]; & 运行结果如下： &
需要说明的是： 还可以通过搭建simulink模块的方法进行求解；搭建好模块后，设置“configuration parameters”相关项目：在command区输入 t_final=100; x0=[0;0;1e-10]&&& 然后就可以运行仿真了，但仿真完后还看不到结果，需在命令口中给出 plot(tout,yout),plot(tout,yout(:,1)),plot(tout,yout(:,2)),plot(tout,yout(:,3)),plot3(yout(:,1),yout(:,2),yout(:,3)); axis([10 40 -20 20 -20 20]); & 特别注意的是：后面4个ode函数，可以解变系数的方程组，当然如果变量很多，系数复杂的话，求解时间会相当长。呵呵
阅读(4807)|
用微信&&“扫一扫”
将文章分享到朋友圈。
用易信&&“扫一扫”
将文章分享到朋友圈。
历史上的今天
loftPermalink:'',
id:'fks_',
blogTitle:'[转载]matlab 微分方程数值解函数',
blogAbstract:'
${a.selfIntro|escape}{if great260}${suplement}{/if}
{list a as x}
推荐过这篇日志的人：
{list a as x}
{if !!b&&b.length>0}
他们还推荐了：
{list b as y}
转载记录：
{list d as x}
{list a as x}
{list a as x}
{list a as x}
{list a as x}
{if x_index>4}{break}{/if}
${fn2(x.publishTime,'yyyy-MM-dd HH:mm:ss')}
{list a as x}
{if !!(blogDetail.preBlogPermalink)}
{if !!(blogDetail.nextBlogPermalink)}
{list a as x}
{if defined('newslist')&&newslist.length>0}
{list newslist as x}
{if x_index>7}{break}{/if}
{list a as x}
{var first_option =}
{list x.voteDetailList as voteToOption}
{if voteToOption==1}
{if first_option==false},{/if}&&“${b[voteToOption_index]}”&&
{if (x.role!="-1") },“我是${c[x.role]}”&&{/if}
&&&&&&&&${fn1(x.voteTime)}
{if x.userName==''}{/if}
网易公司版权所有&&
{list x.l as y}
{if defined('wl')}
{list wl as x}{/list}小木虫 --- 700万学术达人喜爱的学术科研平台
热门搜索：
&&【原创】微分方程求解用什么软件好
【原创】微分方程求解用什么软件好
我想解常微分方程，老板讲VC比较好，一定要用VC编写程序，我看到很多用MATLAB,，我认为MATLAB比较好。不知道学物理和数学的研究生都用什么程序编程？
学术必备与600万学术达人在线互动！
扫描下载送金币
北京学而思教育科技有限公司 地址：北京市海淀区北三环甲18号中鼎大厦A座1层102室 电话：010-}