matlab数值求解边界条件的微分方程组解一个方程组,文献上给出了如下的形式:dx/dz=-a1*x-g1*v1/v2*x(m+n+2d)+c1*ydy/dz=a1*y-g1*v1/v2*x(m+n+2d)-c1*xdm/dz=-a2*m+g1*(m+d)*(x+y)-g2*v2/vs*m*(w+v)+c2*ndn/dz=a2*n-g1*(n+d)*(x+y)+g2*v
来源:学生作业帮助网 编辑:六六作业网 时间:2024/12/23 11:02:28
matlab数值求解边界条件的微分方程组解一个方程组,文献上给出了如下的形式:dx/dz=-a1*x-g1*v1/v2*x(m+n+2d)+c1*ydy/dz=a1*y-g1*v1/v2*x(m+n+2d)-c1*xdm/dz=-a2*m+g1*(m+d)*(x+y)-g2*v2/vs*m*(w+v)+c2*ndn/dz=a2*n-g1*(n+d)*(x+y)+g2*v
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.
matlab数值求解边界条件的微分方程组解一个方程组,文献上给出了如下的形式:dx/dz=-a1*x-g1*v1/v2*x(m+n+2d)+c1*ydy/dz=a1*y-g1*v1/v2*x(m+n+2d)-c1*xdm/dz=-a2*m+g1*(m+d)*(x+y)-g2*v2/vs*m*(w+v)+c2*ndn/dz=a2*n-g1*(n+d)*(x+y)+g2*v
请核实一下条件:
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,根本没法往下算.