Question

matlab数值求解边界条件的微分方程组

解一个方程组，文献上给出了如下的形式：

dx/dz=-a1*x-g1*v1/v2*x(m+n+2d)+c1*y
dy/dz=a1*y-g1*v1/v2*x(m+n+2d)-c1*x

dm/dz=-a2*m+g1*(m+d)*(x+y)-g2*v2/vs*m*(w+v)+c2*n
dn/dz=a2*n-g1*(n+d)*(x+y)+g2*v2/vs*n*(w+v)-c2*m
dw/dz=as*v-g2*v2/vs*v*(m+n)-gb*w*v
dv/dz=-as*w+g2*v2/vs*w*(m+n)-gb*w*v
其中z是自变量，x y m n w v 是函数，其他的都是常数。
边界条件和初值条件为：

x(0)=0.47 y(50)=0.47 w(0)=0.0075 v(0)=0.7e-6
不知道4个条件能不能做。。。
按照matlab求解边界条件微分方程bvp4c方法，我进行了如下处理：
y(1)=x,y(2)=y,y(3)=m,y(4)=n,y(5)=w,y(6)=v
方程M文件
function dydx=ivode(x,y)
dydx=[-a1*y(1)-g1*v1/v2*y(1)*(y(3)+y(4)+2*d)+c1*y(2);
a1*y(2)+g1*v1/v2*y(2)*(y(3)+y(4)+2*d)-c1*y(1);
-a2*y(3)+g1*(y(3)+d)*(y(1)+y(2))-g2*v2/vs*y(3)*(y(5)+y(6))+c2*y(4);
a2*y(4)-g1*(y(4)+d)*(y(1)+y(2))+g2*v2/vs*y(4)*(y(5)+y(6))-c2*y(3);
as*y(6)-g2*v2/vs*y(6)*(y(3)+y(4))-gb*y(5)*y(6);
-as*y(5)+g2*v2/vs*y(5)*(y(3)+y(4))-gb*y(5)*y(6)];
边界条件x(0)=0.47 y(50)=0.47 w(0)=0.0075 v(50)=0.7e-6
边界条件M文件

function res=ivbc(ya,yb)
res=[ya(1)-0.47;ya(5)-0.0075;yb(2)-0.47;yb(6)-0.7e-6];
各种常数参数如下：
a1=0.31;a2=0.25;as=0.2;
g1=0.51;g2=0.38;gb=5e-14;
h=6.62e-34;kb=1.38e-23;t=298;
c1=1.0e-5;c2=6.0e-5;
v1=3e8/1365e-9;v2=3e8/1455e-9;vs=3e8/1550e-9;
dv2=0.4e-12;
d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)*kb*t)-1));

然后再怎么做，如果需要6个边界条件的话，暂定ya(3)=0,yb(4)=0。求高手指点啊，最好有完整的程序，自己调了半天全是错。。。
1楼的大神指出的错误是我自己没打对，修改一下
d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)/(kb*t))-1));

算出d的值为1.2298e-31。

Answer 1

请核实一下条件：
d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)*kb*t)-1));
这个式子确定无误吗？

按照现在给的数值：h=6.62e-34;kb=1.38e-23;
h*(v1-v2)*kb*t的值约为3.7e-041，所以exp(h*(v1-v2)*kb*t)-1=0，d为inf，根本没法往下算。

追问

感谢大神指出错误，是我自己打错了，现在更正过来，
d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)/(kb*t))-1)); 希望大神指点一下接下怎么做，看了不少MATLAB的书，很少有教怎么解边界条件下微分方程的，求指点啊。

追答

function zd513599949
% 解的初始估计
solinit = bvpinit(linspace(0,50,5),[1 0 0 0 0 0]);

% BVP问题求解
sol = bvp4c(@ivode,@ivbc,solinit);

% 结果绘图
t=sol.x;
vars={'x', 'y', 'm', 'n', 'w', 'v'};
for i=1:length(vars)
subplot(2,3,i);
plot(t,sol.y(i,:));
xlabel('t');
ylabel(vars{i});
end

function dydx=ivode(x,y)

% 常数定义
a1=0.31;a2=0.25;as=0.2;
g1=0.51;g2=0.38;gb=5e-14;
h=6.62e-34;kb=1.38e-23;t=298;
c1=1.0e-5;c2=6.0e-5;
v1=3e8/1365e-9;v2=3e8/1455e-9;vs=3e8/1550e-9;
dv2=0.4e-12;
% d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)*kb*t)-1))
d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)/(kb*t))-1));

% 微分方程
dydx=[-a1*y(1)-g1*v1/v2*y(1)*(y(3)+y(4)+2*d)+c1*y(2);
a1*y(2)+g1*v1/v2*y(2)*(y(3)+y(4)+2*d)-c1*y(1);
-a2*y(3)+g1*(y(3)+d)*(y(1)+y(2))-g2*v2/vs*y(3)*(y(5)+y(6))+c2*y(4);
a2*y(4)-g1*(y(4)+d)*(y(1)+y(2))+g2*v2/vs*y(4)*(y(5)+y(6))-c2*y(3);
as*y(6)-g2*v2/vs*y(6)*(y(3)+y(4))-gb*y(5)*y(6);
-as*y(5)+g2*v2/vs*y(5)*(y(3)+y(4))-gb*y(5)*y(6)];

% 边界条件
function res=ivbc(ya,yb)
res=[ya(1)-0.47; ya(5)-0.0075; yb(2)-0.47; ya(6)-0.7e-6; ya(3); yb(4)];