您的位置:首页 > 编程语言 > MATLAB

求解Matlab微分方程组中的时移问题!!!

2009-12-27 19:26 381 查看
问题来源:http://www.ilovematlab.cn/thread-61233-1-1.html

问题描述:

dx/dt=-x*e^(1-t)+0.8y(t-0.1). x(0)=0.
dy/dt=x-y^3. y(0)=2.
上面给定的是微分方程组及其初始条件,要求作出两个图形来。
在这我想做点补充,有些朋友可能会误解我这个方程组,其中0.8y(t-0.1)中的 y(t-0.1)不是y*(t-0.1),而是函数y(t)向右移了0.1个单位即y(t-0.1).
Forcal代码:
!using["XSLSF"];                //使用命名空间XSLSF
//数组xArray存放x的值;ti为当前有效值的个数;tmax为ti对应的时间;tmin为起始时间。
xt(t:k:xArray,ti,tmax,tmin)=
{
k=(t-tmin)/(tmax-tmin)*ti-1,
if[k<0, k=0], if[k>=ti, k=ti-1],
xArray.getrai[k]
};
//数组yArray存放y的值;ti为当前有效值的个数;tmax为ti对应的时间;tmin为起始时间。
yt(t:k:yArray,ti,tmax,tmin)=
{
if{tmax==tmin,return[yArray.getrai(0)]},
k=(t-tmin)/(tmax-tmin)*ti-1,
if[k<0, k=0], if[k>=ti, k=ti-1],
yArray.getrai[k]
};
//函数定义,用于计算微分方程组中各方程右端函数值,连分式法对微分方程组积分一步函数pbs1将调用该函数f。
//t为自变量,x,y为函数值,dx,dy为右端函数值(即微分方程的值)。
f(t,x,y,dx,dy)=
{
dx=-x*exp(1-t)+0.8*yt(t-0.1),
dy=x-y^3
};
//用连分式法对微分方程组进行积分,获得理论值。
//t1,t2为积分的起点和终点。
//h,s为自动变量。
//模块变量:hf为函数f的句柄,要预先获得该句柄;Array为工作数组;step为积分步长;eps为积分精度。
积分(t1,t2:h,s:hf,Array,step,eps)=
{
s=ceil[(t2-t1)/step],        //计算积分步数
h=(t2-t1)/s,                 //重新计算积分步长
{   pbs1[hf,t1,Array,h,eps], //对微分方程组积分一步
t1=t1+h                  //积分步长增加
}.until[abs(t1-t2)<h/2]      //连续积分至t2
};
数据(:i,h,t1,t2,x,y,static,free:hf,Array,step,eps,max,xArray,yArray,ti,tmax,tmin)= //微分方程组积分获得数据
{
if[free,delete(Array),delete(xArray),delete(yArray),return(0)],
hf=HFor("f"),                               //模块变量hf保存函数f的句柄,预先用函数HFor获得该句柄
step=0.01,eps=1e-6,                         //积分步长step要合适,积分精度eps越小越精确,用于对微分方程组积分一步函数pbs1
Array=new[rtoi(real_s),rtoi(2*15)],         //申请工作数组
max=1001,
xArray=new[rtoi(real_s),rtoi(max)],         //申请工作数组
yArray=new[rtoi(real_s),rtoi(max)],         //申请工作数组
Array.setra(0,0,2),                         //设置积分初值,通过模块变量Array传递,Array是一个数组
xArray.setra(0,0),                          //设置xArray的第一个值
yArray.setra(0,2),                          //设置yArray的第一个值
ti=1, h=step*2, tmin=0, tmax=0, t1=0, t2=h,
i=1,(i<max).while{
积分(&t1,t2),                             //从t1积分到t2
Array.getra(0,&x,&y),                     //从数组Array获得t2时的积分值
xArray.setra(i,x), yArray.setra(i,y),     //将积分值保存到数组
ti=i+1, tmax=t1, t2=t2+h,
i++
}
};
//绘制函数图形
(::hxt,hyt)= hxt=HFor("xt"), hyt=HFor("yt");  //获得函数xt和yt的句柄,保存在模块变量中
gleDrawScene[HFor("Scene")],stop();      //设置场景绘制函数后退出
Scene(::hxt,hyt,tmax,tmin)=
{
glClear[],                           //清除屏幕以及深度缓存
glLoadIdentity[],                    //重置视图
glTranslated[0,0,-20],               //移动坐标,向屏幕里移动10个单元
glColor3d[0,0,1],                    //设置颜色
fgPlot[hxt,tmin,tmax,FG_Y,0,10],     //绘制一元函数图象
glColor3d[1,0,0],                    //设置颜色
fgPlot[hyt,tmin,tmax,FG_Y,0,10]      //绘制一元函数图象
};
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: