|
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ 感谢参与,应助指数 +1 亮亮1283: 金币+20, ★★★★★最佳答案, 哈哈,就是这个原因。太谢谢你了! 2012-05-28 15:17:02
不知道对不对,画出来的结果好像都是1*10^16
原程序主要有两个问题,你的那些参数需要用全局变量,子函数不能直接用主函数的值的
还有xlabel,ylabel,zlabel,不是xlable等
function pdex1()
clc
clear
global Up Kb T P a y r h v Un e tp p0 n0 n1 p1 G tn
Up =400;
Un =8000;
Kb =1.38*10^(-23);
h =6.63*10^(-34);
e =1.6*10^(-19);
r =0.35;
y =0.6;
n0 =1*10^16;
p0 =1*10^16;
tn =1*10^(-9);
tp =1*10^(-9);
p1 =p0;
n1 =n0;
T =300;
a =1.65*10^4;
v =(3/532)*10^17;
X =10;
P=10;
%主函数----------------------------------------------------------------
m=0;
x=0:0.1:10;
t=0:20:1000;
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% u =[ sol(:,:,1), sol(:,:,2)];
figure('numbertitle','off','name','PDE Demo-by matlab')
subplot(211)
surf(x,t,sol(:,:,1))
title('The Solution of p')
xlabel('x')
ylabel('t')
zlabel('p')
subplot(212)
surf(x,t,sol(:,:,2))
title('The Solution of n')
xlabel('x')
ylabel('t')
zlabel('n')
%目标函数----------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
global Up Kb T P a y r h v Un e tp p0 n0 n1 p1 G tn
c= [1;1];
f= [(Up*Kb*T/e)*DuDx(1); (Un*Kb*T/e)*DuDx(2)];
G= a*y*(1-r)*P*exp(-a*x)/(h*v);
R1= u(1)*u(2)-p0*n0;
R2=tn*(u(1)+p1)+tp*(u(2)+n1);
R=R1/R2;
s= [1;1].*(G-R);
%初值条件函数------------------------------------------------------------
function u0 = pdex1ic(x)
global p0 n0
u0= [p0;n0];
%边界条件函数------------------------------------------------------------
function [pa,qa,pb,qb] = pdex1bc(xa,ua,xb,ub,t) %a为下边界, b为上边界
pa =[0; 0];
qa =[1; 1];
pb =[0; 0];
qb =[1; 1]; |
|