关于微分方程组的求解！

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';
[N2,pp1,pp2,ps1,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=[0:0.01:16];
%定义各常量的数值
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';
[N2,pp1,pp2,ps1,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([var_set ' intersect ' var_set]);

Error in ==> m_10_24 at 29
[N2,pp1,pp2,ps1,ps2]=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=[0:0.01:16];
%定义各常量的数值
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

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

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

不好解，因为第一个方程是个代数方程。这是一个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[hf,t1,Array,h,eps], //对微分方程组积分一步
t1=t1+h //积分步长增加
}.until[abs(t1-t2)<h/2] //连续积分至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[rtoi(real_s),rtoi(4*15)],
x=new[rtoi(real_s),rtoi(4),rtoi(EndType),1,1],
i=netn[HFor("求初值"),1e-15,0.1,0.1,x,100],
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[Array],delete[x]
};
验证(:pp20,ps20,t1,t2,s1,s2,s3,s4,h:hf,Array,step,eps,r1)=
{
Array=new[rtoi(real_s),rtoi(4*15)],
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[hf,t1,Array,h,eps], //对微分方程组积分一步
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[abs(t1-t2)<h/2], //连续积分至t2
delete[Array]
};

复制代码

结果：

迭代次数=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

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

账号		自动登录	找回密码
密码			注册

账号		自动登录	找回密码
密码			立即注册

关于微分方程组的求解！

本帖子中包含更多资源

评分

相关帖子

关于常微分方程的求解？

评分

本帖子中包含更多资源

评分