论坛 › K01：MatLab › 关于微分方程组的求解！

ln_buaa 发表于 2009-10-25 15:55:15

关于微分方程组的求解！

这个方程组，最初用的dsolve命令，总报错！源程序如下：
clear
xap=6*exp(-25);
xep=2.5*exp(-26);
tp=0.00474;
ts=0.82;
h=6.626*exp(-34);
vp=3.07*exp(14);
A=pi*(100*exp(-12));
vs=2.76*exp(14);
xas=1.4*exp(-27);
xes=2*exp(-25);
t1=1*exp(-3);
N=0.9*exp(25);
ap=3*exp(-3);
as=5*exp(-3);
r1=0.998;
r2=0.04;
%定义变量
syms N2 pp1 pp2 ps1 ps2 t
%定义方程组中各等式
eq1='N2=N*((pp1*xap*tp/(h*vp*A)+(ts/(h*vs*A))*xas*(ps1+ps2))/((pp1*(xap+xep)*tp)/h*vp*A+1/t1+ts/(h*vs*A)*(xes+xas)*(ps1+ps2)))';
eq2='D(pp1)=-tp*(xap*N-(xap+xep)*N2)*pp1-ap*pp1';
eq3='-D(pp2)=-tp*(xap*N-(xap+xep)*N2)*pp2-ap*pp2';
eq4='D(ps1)=ts*((xes+xas)*N2-xas*N)*ps1-as*ps1';
eq5='-D(ps2)=ts*((xes+xas)*N2-xas*N)*ps2-as*ps2';
=dsolve('eq1','eq2','eq3','eq4','eq5','pp1(0)=2000','pp2(16)=0','ps1(0)=r1*ps2(0)','ps2(16)=r2*ps1(16)')
但是求不到方程的封闭解，如果用ode45命令求解，如何改写程序，最近刚学MATLAB，请牛牛们指教下，不胜感激

ln_buaa 发表于 2009-10-25 10:10:21

关于常微分方程的求解？

lear
%定义变量z的取值区间
t=;
%定义各常量的数值
xap=6*exp(-25);
xep=2.5*exp(-26);
tp=0.00474;
ts=0.82;
h=6.626*exp(-34);
vp=3.07*exp(14);
A=pi*(100*exp(-12));
vs=2.76*exp(14);
xas=1.4*exp(-27);
xes=2*exp(-25);
t1=1*exp(-3);
N=0.9*exp(25);
ap=3*exp(-3);
as=5*exp(-3);
r1=0.998;
r2=0.04;
%定义变量
syms N2 pp1 pp2 ps1 ps2 t
%定义方程组中各等式
eq1='N2/N=(pp1*xap*tp/(h*vp*A)+(ts/(h*vs*A))*xas*(ps1+ps2))/((pp1*(xap+xep)*tp)/h*vp*A+1/t1+ts/(h*vs*A)*(xes+xas)(ps1+ps2))';
eq2='D(pp1)=-tp*(xap*N-(xap+xep)*N2)*pp1-ap*pp1';
eq3='-D(pp2)=-tp*[xap*N-(xap+xep)*N2)*pp2-ap*pp2';
eq4='D(ps1)=ts*((xes+xas)*N2-xas*N)*ps1-as*ps1';
eq5='-D(ps2)=ts*((xes+xas)*N2-xas*N)*ps2-as*ps2';
=dsolve('eq1','eq2','eq3','eq4','eq5','pp1(0)=2000','pp2(16)=0','ps1(0)=r1*ps2(0)','ps2(16)=r2*ps1(16)')
这个不一定能接出来，但是出现了这样的报错：
??? Error using ==> maple at 129
at offset 2, `}` unexpected

Error in ==> dsolve at 141
var_set = maple();

Error in ==> m_10_24 at 29
=dsolve('eq1','eq2','eq3','eq4','eq5','pp1(0)=2000','pp2(16)=0','ps1(0)=r1*ps2(0)','p
因为开始接触MATLAB的时间不是很长，不知道错在哪了？
哪位经验丰富的牛牛，指点下

TBE_Legend 发表于 2009-10-25 14:02:15

lear
%定义变量z的取值区间
t=;
%定义各常量的数值
xap=6*exp(-25);
xep=2.5*exp(-26);
tp=0.00474;
ts=0.82;
h=6.626*exp(-34);
vp=3.07*exp(14);
A=pi*(100*exp(-12));
vs=2.76*exp(14);
x ...
ln_buaa 发表于 2009-10-25 10:10 http://forum.simwe.com/images/common/back.gif

能不能采用数学描述把你的问题说下？ 其中的一些公式写到文本文件中，以附件的形式穿上来。

TBE_Legend 发表于 2009-10-25 18:27:15

这个方程组，最初用的dsolve命令，总报错！源程序如下：
clear
xap=6*exp(-25);
xep=2.5*exp(-26);
tp=0.00474;
ts=0.82;
h=6.626*exp(-34);
vp=3.07*exp(14);
A=pi*(100*exp(-12));
vs=2.76*exp(14);
xas ...
ln_buaa 发表于 2009-10-25 15:55 http://forum.simwe.com/images/common/back.gif

不好解，因为第一个方程是个代数方程。这是一个DAE。

1）能否把你的第一个方程改变下，使之出现N2的导数项。
2）要解DAE的话，一般是初值问题比较好办，所以如果1）不可行的话，你得把你的边界条件改改，使之全部为初始条件。

wanglu 发表于 2009-10-25 21:47:12

用Forcal试了一下，先求两个未知初值，再求其他数据。
如果公式没有错误，感觉数据不是很理想：


//这里是代码窗口，请将Forcal代码写在下面
!using["XSLSF"]; //使用命名空间XSLSF
!(::pi,xap,xep,tp,ts,h,vp,A,vs,xas,xes,t1,N,ap,as,r1,r2)= //定义常量
{
pi=acos(0)*2,
xap=6*exp(-25),
xep=2.5*exp(-26),
tp=0.00474,
ts=0.82,
h=6.626*exp(-34),
vp=3.07*exp(14),
A=pi*(100*exp(-12)),
vs=2.76*exp(14),
xas=1.4*exp(-27),
xes=2*exp(-25),
t1=1*exp(-3),
N=0.9*exp(25),
ap=3*exp(-3),
as=5*exp(-3),
r1=0.998,
r2=0.04
};
//定义方程
N2(pp1,pp2,ps1,ps2::xap,xep,tp,ts,h,vp,A,vs,xas,xes,t1,N,ap,as,r1,r2)=N*((pp1*xap*tp/(h*vp*A)+(ts/(h*vs*A))*xas*(ps1+ps2))/((pp1*(xap+xep)*tp)/h*vp*A+1/t1+ts/(h*vs*A)*(xes+xas)*(ps1+ps2)));
f(t,pp1,pp2,ps1,ps2,Dpp1,Dpp2,Dps1,Dps2::tp,xap,xep,N,ap,ts,xes,xas,as)=
{
Dpp1=-tp*(xap*N-(xap+xep)*N2(pp1,pp2,ps1,ps2))*pp1-ap*pp1,
Dpp2=-{-tp*(xap*N-(xap+xep)*N2(pp1,pp2,ps1,ps2))*pp2-ap*pp2},
Dps1=ts*((xes+xas)*N2(pp1,pp2,ps1,ps2)-xas*N)*ps1-as*ps1,
Dps2=-{ts*((xes+xas)*N2(pp1,pp2,ps1,ps2)-xas*N)*ps2-as*ps2}
};
getL(t1,t2:h:hf,Array,step,eps)=
{
    h=(t2-t1)/step,            //计算积分步长
    {   pbs1, //对微分方程组积分一步
      t1=t1+h                  //积分步长增加
    }.until      //连续积分至t2
};
求初值(pp20,ps20,y1,y2:t1,t2,t3,t4:Array,r1,r2)=    //函数定义
{
Array.setra(0,2000,pp20,r1*ps20,ps20),//设置积分初值，通过模块变量Array传递，Array是一个数组
getL(0,16),                           //从0积分到16
Array.getra(0,&t1,&t2,&t3,&t4),         //从数组Array获得t2时的积分值
y1=t2,
y2=t4-r2*t3
};
main(:x,t1,t2,i:hf,Array,step,eps)=
{
hf=HFor("f"),step=50,eps=1e-9,
Array=new,
x=new,
i=netn,
x.getra(0,&t1,&t2),
//若返回值（迭代次数）为-1表示代数方程奇异，若返回值为-2表示β=0，这两种情况可放宽精度要求、改变各个初值或改变各个方程顺序再试；
//若返回值等于0说明迭代了k次仍未满足精度要求，程序工作失败。
printff{"迭代次数={1,i}, pp20={2,r}, ps20={3,r}\r\n",i,t1,t2},
delete,delete
};
验证(:pp20,ps20,t1,t2,s1,s2,s3,s4,h:hf,Array,step,eps,r1)=
{
    Array=new,
    pp20=-3.5296322379324667e-023, ps20=1.8589081085716699e-019,
    Array.setra(0,2000,pp20,r1*ps20,ps20),//设置积分初值，通过模块变量Array传递，Array是一个数组
    t1=0,t2=16,step=50,eps=1e-9,
    h=(t2-t1)/step,            //计算积分步长
    {   pbs1, //对微分方程组积分一步
      t1=t1+h,               //积分步长增加
      Array.getra(0,&s1,&s2,&s3,&s4),         //从数组Array获得t2时的积分值
      printff{"t={1,r,-10.5}n2={2,r,-18.8}pp1={3,r,-18.8}pp2={4,r,-18.8}ps1={5,r,-18.8}ps2={6,r,-18.8}\r\n",t1,N2(s1,s2,s3,s4),s1,s2,s3,s4}
    }.until,   //连续积分至t2
    delete
};


结果：

迭代次数=3, pp20=-3.5296322379324667e-023, ps20=1.8589081085716699e-019
0.
t=0.32      n2=1106.4052         pp1=1891.1042         pp2=-3.7328797e-023   ps1=1.6381756e-019    ps2=2.1051642e-019   
t=0.64      n2=1106.4051         pp1=1788.1375         pp2=-3.9478309e-023   ps1=1.4465466e-019    ps2=2.3840426e-019   
t=0.96      n2=1106.4051         pp1=1690.7771         pp2=-4.1751596e-023   ps1=1.2773338e-019    ps2=2.699865e-019   
t=1.28      n2=1106.4051         pp1=1598.7178         pp2=-4.4155786e-023   ps1=1.127915e-019   ps2=3.0575256e-019   
t=1.6       n2=1106.4051         pp1=1511.671          pp2=-4.6698418e-023   ps1=9.9597481e-020    ps2=3.4625667e-019   
t=1.92      n2=1106.4051         pp1=1429.3637         pp2=-4.9387462e-023   ps1=8.7946859e-020    ps2=3.9212649e-019   
t=2.24      n2=1106.4051         pp1=1351.5378         pp2=-5.223135e-023    ps1=7.7659092e-020    ps2=4.4407286e-019   
t=2.56      n2=1106.4051         pp1=1277.9494         pp2=-5.5238998e-023   ps1=6.8574759e-020    ps2=5.0290074e-019   
t=2.88      n2=1106.4051         pp1=1208.3677         pp2=-5.8419836e-023   ps1=6.0553084e-020    ps2=5.6952176e-019   
t=3.2       n2=1106.4051         pp1=1142.5746         pp2=-6.1783837e-023   ps1=5.3469761e-020    ps2=6.4496829e-019   
t=3.52      n2=1106.4051         pp1=1080.3638         pp2=-6.5341547e-023   ps1=4.7215025e-020    ps2=7.3040949e-019   
t=3.84      n2=1106.4051         pp1=1021.5402         pp2=-6.9104123e-023   ps1=4.1691949e-020    ps2=8.2716939e-019   
t=4.16      n2=1106.4051         pp1=965.91949         pp2=-7.3083359e-023   ps1=3.6814947e-020    ps2=9.367474e-019   
t=4.48      n2=1106.4051         pp1=913.32719         pp2=-7.7291732e-023   ps1=3.2508442e-020    ps2=1.0608416e-018   
t=4.8       n2=1106.4051         pp1=863.59843         pp2=-8.1742437e-023   ps1=2.87057e-020      ps2=1.201375e-018   
t=5.12      n2=1106.4051         pp1=816.57729         pp2=-8.6449428e-023   ps1=2.5347791e-020    ps2=1.3605253e-018   
t=5.44      n2=1106.4051         pp1=772.11636         pp2=-9.1427462e-023   ps1=2.238268e-020   ps2=1.5407589e-018   
t=5.76      n2=1106.4051         pp1=730.07624         pp2=-9.6692148e-023   ps1=1.976442e-020   ps2=1.7448686e-018   
t=6.08      n2=1106.4051         pp1=690.32511         pp2=-1.0225999e-022   ps1=1.7452436e-020    ps2=1.9760174e-018   
t=6.4       n2=1106.4051         pp1=652.73834         pp2=-1.0814845e-022   ps1=1.5410901e-020    ps2=2.2377873e-018   
t=6.72      n2=1106.4051         pp1=617.1981          pp2=-1.1437598e-022   ps1=1.3608179e-020    ps2=2.5342348e-018   
t=7.04      n2=1106.4051         pp1=583.59295         pp2=-1.2096212e-022   ps1=1.2016334e-020    ps2=2.8699537e-018   
t=7.36      n2=1106.4051         pp1=551.81753         pp2=-1.279275e-022    ps1=1.0610698e-020    ps2=3.2501465e-018   
t=7.68      n2=1106.405          pp1=521.77222         pp2=-1.3529397e-022   ps1=9.3694897e-021    ps2=3.6807048e-018   
t=8.      n2=1106.405          pp1=493.36281         pp2=-1.4308463e-022   ps1=8.2734742e-021    ps2=4.1683006e-018   
t=8.32      n2=1106.405          pp1=466.50023         pp2=-1.513239e-022    ps1=7.3056674e-021    ps2=4.7204899e-018   
t=8.64      n2=1106.405          pp1=441.10027         pp2=-1.6003762e-022   ps1=6.4510718e-021    ps2=5.3458297e-018   
t=8.96      n2=1106.405          pp1=417.08328         pp2=-1.692531e-022    ps1=5.6964442e-021    ps2=6.0540104e-018   
t=9.28      n2=1106.405          pp1=394.37396         pp2=-1.7899923e-022   ps1=5.0300908e-021    ps2=6.8560063e-018   
t=9.6       n2=1106.405          pp1=372.90112         pp2=-1.8930658e-022   ps1=4.4416855e-021    ps2=7.7642453e-018   
t=9.92      n2=1106.405          pp1=352.59743         pp2=-2.0020746e-022   ps1=3.9221101e-021    ps2=8.792802e-018   
t=10.24   n2=1106.405          pp1=333.39924         pp2=-2.1173604e-022   ps1=3.4633132e-021    ps2=9.9576152e-018   
t=10.56   n2=1106.4049         pp1=315.24634         pp2=-2.2392848e-022   ps1=3.0581849e-021    ps2=1.1276735e-017   
t=10.88   n2=1106.4049         pp1=298.08184         pp2=-2.36823e-022   ps1=2.7004474e-021    ps2=1.2770604e-017   
t=11.2      n2=1106.4049         pp1=281.8519          pp2=-2.5046002e-022   ps1=2.3845569e-021    ps2=1.446237e-017   
t=11.52   n2=1106.4049         pp1=266.50565         pp2=-2.6488231e-022   ps1=2.1056184e-021    ps2=1.6378251e-017   
t=11.84   n2=1106.4049         pp1=251.99498         pp2=-2.8013508e-022   ps1=1.8593094e-021    ps2=1.8547935e-017   
t=12.16   n2=1106.4049         pp1=238.27438         pp2=-2.9626615e-022   ps1=1.6418128e-021    ps2=2.1005044e-017   
t=12.48   n2=1106.4048         pp1=225.30083         pp2=-3.1332611e-022   ps1=1.4497583e-021    ps2=2.3787656e-017   
t=12.8      n2=1106.4048         pp1=213.03367         pp2=-3.3136842e-022   ps1=1.2801698e-021    ps2=2.693889e-017   
t=13.12   n2=1106.4048         pp1=201.43443         pp2=-3.5044968e-022   ps1=1.1304193e-021    ps2=3.0507578e-017   
t=13.44   n2=1106.4048         pp1=190.46675         pp2=-3.7062969e-022   ps1=9.9818616e-022    ps2=3.4549024e-017   
t=13.76   n2=1106.4048         pp1=180.09623         pp2=-3.9197173e-022   ps1=8.8142126e-022    ps2=3.9125854e-017   
t=14.08   n2=1106.4047         pp1=170.29037         pp2=-4.1454272e-022   ps1=7.7831518e-022    ps2=4.4308992e-017   
t=14.4      n2=1106.4047         pp1=161.01841         pp2=-4.3841341e-022   ps1=6.8727015e-022    ps2=5.0178759e-017   
t=14.72   n2=1106.4047         pp1=152.25129         pp2=-4.6365866e-022   ps1=6.0687529e-022    ps2=5.6826114e-017   
t=15.04   n2=1106.4047         pp1=143.96153         pp2=-4.903576e-022    ps1=5.358848e-022   ps2=6.4354068e-017   
t=15.36   n2=1106.4046         pp1=136.12312         pp2=-5.1859396e-022   ps1=4.7319856e-022    ps2=7.2879276e-017   
t=15.68   n2=1106.4046         pp1=128.7115          pp2=-5.4845626e-022   ps1=4.1784517e-022    ps2=8.2533848e-017   
t=16.       n2=1106.4046         pp1=121.70343         pp2=-5.8003812e-022   ps1=3.6896687e-022    ps2=9.3467394e-017

wanglu 发表于 2009-10-25 21:59:39

感觉数据不是很理想：试了多次，两个未知初值总在0附近？

TBE_Legend 发表于 2009-10-25 22:39:18

本帖最后由 TBE_Legend 于 2009-10-26 00:48 编辑

感觉数据不是很理想：试了多次，两个未知初值总在0附近？
wanglu 发表于 2009-10-25 21:59 http://forum.simwe.com/images/common/back.gif

pp1红色，怎么其他的都是0附近？难道我数据抄错了。

查看完整版本: 关于微分方程组的求解！