matlab pdebc,matlab的pdepe函数求解，其结果为何只有一个初始时间的解

本帖最后由 ph0099 于 2020-12-17 10:13 编辑

各位大佬，我用pdepe函数求解微分方程，为何在第一次求解之后，其后所有的结果都只有初始时刻的值。

也就是说，我第一次用pdepe函数求解方程得到一个11*120维的解sol，然后我把最后一个时间点的结果sol(11,:)作为初值再求解这个方程，得到的sol就只有1*120，也就是只有初始时刻的值，其后面时刻的结果就没有了，这是为什么？

代码如下

N1 = test_n_new_1(0.0989,80e-6,10e-6);

function N = test_n_new_1(VM,During,Risetime)

% warning off

global Vm AP sigma Ni1 i

vm=zeros(2,1);

vm(2,1)=VM;

% load("n","Ni")

% N(1,1) =  1.5e9;

x=linspace(0.65e-9,15e-9,120);

Ni1 = zeros(1,size(x,2));

AP = zeros(1,1);

N = zeros(1,size(x,2));

step = 2e-6;

gen = 2;

for T=0:step:step %During + 2 * Risetime

Vm=vm(gen-1,1);%

%     Vm=VM;%

sigma1= 2e-2;

sigma0= 1e-3;

sigma = 2*sigma1-(2*sigma1-sigma0)/(1-AP(1,1));

%     sigma =9.9856e-6;

m=0;

Ni1 = N(gen-1,:);

%     Ni1 = Ni;

tfspan = linspace(T,T+step,3);

i=1;

sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,tfspan);

N(gen,:) = sol(end,:);

AP(1,1) = sum(x.^2*pi.*N(gen,:),2);

gen=gen+1;

end

end

function [c,f,s] = pdefun(x,t,u,ux)

global Vm sigma

D=2e-13;

T=310;

k=1.38e-23;    %玻尔兹曼常数

gama=2e-11;

beta=1.4e-19;

rh=0.95e-9;

rt=0.23e-9;

r1=0.65e-9;

Fmax=0.69e-9;

% F = Fmax*Vm^2*((x+rt)/(x + rh + rt));

% E = beta*r1^4/x^4 + 2*pi*gama*x - pi*sigma*x^2 - Fmax*Vm^2*(x + log(x + rh + rt)*(-rh) - log(rh + rt)*(-rh));%integral(fun,0,x)

% dE = -4*beta*r1^4/x^5 + 2*pi*gama -2*pi*x*sigma - Fmax*Vm^2*((x+rt)/(x + rh + rt)) ;

% d2E = (20*beta*r1^4)/x^6 - 2*pi*sigma - (Fmax*Vm^2*rh)/(rh + rt + x)^2;

c=1;

f = D*ux+D*u*(-4*beta*r1^4/x^5 + 2*pi*gama -2*pi*x*sigma - Fmax*Vm^2*((x+rt)/(x + rh + rt)))/k/T;

s = 0;

end

function u0=pdeic(x)   %建立偏微分方程的初始条件函数

global Ni1

% u0= Ni1(1,round((x-0.64e-9)*1e11));

global i

u0= Ni1(1,i);

i=i+1;

% u0 = 0;

end

function [pl,ql,pr,qr]=pdebc(xl,ul,xr,ur,t)   %建立偏微分方程的边界条件函数

global Vm sigma

D=2e-13;

T=310;

k=1.38e-23;    %玻尔兹曼常数

gama=2e-11;

A=1e9;

B=20;%*k*T

v=0.33e3;

Fmax=0.69e-9;

beta=1.4e-19;

rh=0.95e-9;

rt=0.23e-9;

r1=0.65e-9;

pl=A*exp(B*Vm^2/k/T)+D*(-4*beta*r1^4/xl^5 + 2*pi*gama -2*pi*xl*sigma - Fmax*Vm^2*((xl+rt)/(xl+ rh + rt)))*ul/k/T - v*ul;

ql=D;

% ql=D;

pr=D*ur*(-4*beta*r1^4/xr^5 + 2*pi*gama -2*pi*xr*sigma - Fmax*Vm^2*((xr+rt)/(xr + rh + rt)))/k/T;

qr=D;

% qr=D;

end