MATLAB Answers

Problem with solving the inhomogeneous systems of ODE

6 views (last 30 days)
raha ahmadi
raha ahmadi on 31 Jul 2020
Edited: raha ahmadi on 1 Aug 2020
Hi, I 'm trying to solve a system of ODE neumerically, in these equaions I have a funcion of time ,Pc, which was calculated seperately . I dont know how enter this function that now is a matrix of numbers to this system of odes. I know it s better to enter it with time dependent form but its form is complex and I m not able to extract its coeffisient.In fact so far I solve systems with constant inhomogeneous part .
Sorry I have a lot of constants which makes code long
I am new to MATLAB and would really appreciate the help.
With best wishes
Raha
function [t,y] = solveODEPcNew
%% constant
t=linspace(1e-11,-3,1000);
e=1.6*1e-19;
mili=1e-3;
nA=3.2;
nB=3.18;
nP=3.19;
c=3*1e8;
sigma=0.75;
%% structure
R1=24;%k/w
R2=7;
R3=9.4;
micro=1e-6;
lA=300*micro;
wA=micro;
tA=0.15*micro;%tickness
vA=lA*wA*tA;%volume
aA=lA*wA;%area
%pc structure
lP=250*micro;
wP=1.5*micro;
tP=0.25*micro;%tickness
vP=lP*wP*tP;
aP=lP*wP;
%Bragg section
lB=300*micro;
leffB=135*micro;
wB=1.5*micro;
tB=0.25*micro;%tickness
vB=lB*wB*tB;
aB=lB*wB;
RHOv=4.825e3;
V=[vA,vP,vB,vA,vP,vB];
Cv=3.124e2;
Cc=2.46e-7;
Cs=5.32e-5;
Cj=Cv*RHOv.*V;
%average Heat Capacity
C=mean(Cj);
T0=300;
% Ia=70*mili;
% Iast=70*mili;
% IB=20*mili;
% IBst=20*mili;
% Ip=2.7*mili;
% Ipst=2.7*mili;
%% ohmic heat term in sections Pi=a1i*Ii+a2i*Ii^2
rhoN=[9.24e-5,9.24e-5,2.19e-3,2.229e-3];
%thicness of differenr layer of 6 section
dna=[0.15,0.15,2.5,0.5];
dna=dna.*micro;
Ea=0.83;
% a1a=sigma*Ea;
% a2a=sigma^2*sum(rhoN.*dna)/aA;
% Rna=sum(rhoN.*dna)/aA;
%
dnB=[0.25,0.15,2.5,0.5];
dnB=micro.*dnB;
EB=0.95;
a1B=sigma*EB;
a2B=sigma^2*sum(rhoN.*dnB)/aB;
RnB=sum(rhoN.*dnB)/aB;
%
dnp=[0.25,0.15,2.5,0.5];
dnp=micro.*dnp;
Ep=0.95;%Ea=Ea*
l0=lA+lB+lP;%chip &substrate
rho0=2.229e-3;
d0=100*micro;
w0=400*micro;
rc=(rho0*d0)/(w0*l0);
a1a=sigma*Ea;
a2a=sigma^2*sum(rhoN.*dna)/aA;
IB=20*mili;
JB=sigma*IB/(e*vB);
Ba=-5.97*1e-27;
Ntr=1.37*1e24;
Na0=Ntr
NB0=Ntr;
Np0=Ntr;
Naini =3.1642e+24% with 70 mili
NBini=3e24;
Npini=4e24;
nAini=nA+Ba*(Naini-Na0);
nBini=nB+Ba*(NBini-NB0);
nPini=nP+Ba*(Npini-Np0);
Leffcini=nAini*lA+nBini*leffB+nPini*lP;
mNum= 2*Leffcini/(1553*1e-9);
LambdaB=1553*1e-9/(2*nBini);
g = 1.454401546848779e-13;
Sssini =1.7418e+22;
gammaN =3.8180e+08;
dalphadN=1.8*1e-21;%m^2
AprimA=zetaA*dalphadN;%A'a
k1=mNum*LambdaB/(Leffcini*g)*(-gammaPini*Ba+c*AprimA);
k2=mNum*LambdaB*lA/(Leffcini*lP)*(-gammaPini*Ba+c*AprimA)/g+(leffB-mNum*LambdaB)/lP;
Naini =3.1642e+24% with 70 mili
Nass = 3.8658e+24
tauB=1.0921e-08;% linear assumption
Tw=302;
b1=0.28*1e8;
b2=(0.28*1e-16-0.48*1e-17);
b3=(1+(1.6*1e-2)*(Tw-T0))*(0.52*1e-41);
%Rp(Np)
p1=b1;
p2=b2;
p3=b3;
b3=(1+(1.6*1e-2)*(Tw-T0))*(0.52*1e-41);
JBini=b1*NBini+b2*NBini^2+b3*NBini^3;
Jpini=p1*Npini+p2*Npini^2+p3*Npini^3;
Iaini=70*mili;
%Jaini=sigma/e*vA*Iaini
LambdaB=1553*1e-9/(2*nBini);
%Na=Naini+k1*tauB*(JB-JBini)*(1-exp(-t./tauB));
%Ja=Ndota+gammaN.*Na+g.*(Na-Ntr).*Sssini;
Jp=Jpini+k2*(JB-JBini);
k1=mNum*LambdaB/(Leffcini*g)*(-gammaPini*Ba+c*AprimA);
%k2=mNum*LambdaB*lA/(Leffcini*lP)*(-gammaPini*Ba+c*AprimA)/g+(leffB-mNum*LambdaB)/lP;
Na=Naini+k1*tauB*(JB-JBini)*(1-exp(-t./tauB));
Ndota=k1*(JB-JBini)*exp(-t./tauB);
Ja=Ndota+gammaN.*Na+g.*(Na-Na0).*Sssini;
Ia=e*vP.*Ja/sigma;
%
Iamax=0.212;
Ipmax=Ip;
IBmax=IB;
PaBar=a1a*Iamax+a2a*Iamax^2;
Pa=a1a.*Ia+a2a.*Ia.^2;
Iast=(-a1a+(a1a^2+sqrt(a1a^2-4*a2a.*(Pa-2*PaBar))))/(2*a1a);
% Past=a1a*Iast+a2a*Iast^2;
IBst=20*mili;
Ip=2.7*mili;
Ipst=2.7*mili;
I=Ia+Iast+IB+IBst+Ip+Ipst;
%P=Pa+Past+PB+PBst+Pp+Ppst;
Ileak=(1-sigma).*I.^2;
Pc=rc*Ileak.^2;
%% Matrix of coeffeisients
m11=-1/(R2*Cc);
m12=-m11;
m13=1/(Cc*R1);
% m14=0;
m21=1/(R2*Cs);
m22=-(1/R2+1/R3)*(1/Cs);
m23=0;
% m24=0;
m31=6/(R2*Cc);
m32=-6/(R2*Cc);
m33=-6/(Cc*R1)-1/(R1*C);
%% inhomogeneous part
U1=Pc/Cc;
U2=0;
U3=-6*Pc/Cc;%
%% solve the system T0
tspan=[1e-11 1e-3];
y0=[0 0 0 ];
[t,y] = ode15s(@odefcn,tspan,y0);
%% ODE function
function dydt = odefcn(~,y)
dydt=zeros(3,1);
dydt(1)= m11*y(1)+ m12*y(2)+ m13*y(3)+U1;
dydt(2)= m21*y(1)+ m22*y(2)+ m23*y(3)+U2;
dydt(3)= m31*y(1)+ m32*y(2)+ m33*y(3)+U3;
end
end
  3 Comments
raha ahmadi
raha ahmadi on 1 Aug 2020
Dear Walter
Thanks for your help. I realized my problem and tried to change it . I did a long calculation and replace U1 with its equation (in fact I replace U1 and U3 with their equations as a function of time.after that I introduced time ,t, as a 4th input as you suggested.
I really appreciated if you check the new code. Also now I face the other error
"Warning: Failure at t=1.000000e-07. Unable to meet integration tolerances without
reducing the step size below the smallest value allowed (2.117582e-22) at time t.
> In ode15s (line 730)
In solveODEPcNew (line 133) " this means ode15s cant solve it.
Thanks in advance
Raha
sorry I fotgot send the other parts of my code. my main problem is what Walter Roberson said. thanks for time
function [t,y] = solveODEPcNew
%% constant
mili=1e-3;
tauB=1.0921e-08;
AprimA =5.850000000000000e-22;
sigma=0.75;
R1=24;%k/w
R2=7;
R3=9.4;
micro=1e-6;
lA=300*micro;
wA=micro;
tA=0.15*micro;%tickness
vA=lA*wA*tA;%volume
aA=lA*wA;%area
lP=250*micro;
wP=1.5*micro;
tP=0.25*micro;%tickness
vP=lP*wP*tP;
aP=lP*wP;
lB=300*micro;
wB=1.5*micro;
tB=0.25*micro;%tickness
vB=lB*wB*tB;
aB=lB*wB;
V=[vA,vP,vB,vA,vP,vB];
Cv=3.124e2;
Cc=2.46e-7;
Cs=5.32e-5;
RHOv=4.825e3;
Cj=Cv*RHOv.*V;
C=mean(Cj);
Ntr=1.37*1e24;
Na0=Ntr
l0=lA+lB+lP;%chip &substrate
rho0=2.229e-3;
d0=100*micro;
w0=400*micro;
rc=(rho0*d0)/(w0*l0);
x=rc*(1-sigma)^2
Ba=-5.97*1e-27;
c=3e8;
IB=20*mili;
IBst=20*mili;
Ip=2.7*mili;
Ipst=2.7*mili;
Ires=IB+IBst+Ip+Ipst;
Ea=0.83;
a1a=sigma*Ea;
rhoN=[9.24e-5,9.24e-5,2.19e-3,2.229e-3];
dna=[0.15,0.15,2.5,0.5];
dna=dna.*micro;
a2a=sigma^2*sum(rhoN.*dna)/aA;
lambda0=1553e-9
dndomega=1e-16;
nA=3.2;
freq=c/(nA*lambda0);
omega=2*pi*freq;
zetaA=0.325;
zetaB=0.488;
zetaP=0.488;
ngA=nA+dndomega*omega;
dgdNa=4.83*1e-21;
VgA=c/ngA
g=zetaA*dgdNa*VgA
Sssini = 1.648956271484272e+22;
gammaN =3.8180e+08;
mNum =2.805341422150676e+03
Naini =3.1642e+24
tauB=1.0921e-08;
JB =8.333333333333332e+32;
%JBini=b1*NBini+b2*NBini^2+b3*NBini^3;
JBini =4.376928000000000e+32;
DeltaJB=JB-JBini;
Leffcini = 0.002178347614300;
gammaPini =3.339100000000000e+11;
NBini=3e24;
nB=3.18;
Npini=4e24;
NB0=Ntr;
nBini=nB+Ba*(NBini-NB0);
LambdaB=1553*1e-9/(2*nBini);
k1=mNum*LambdaB/(Leffcini*g)*(-gammaPini*Ba+c*AprimA);
T0=300;
Iamax=0.212;
Ipmax=Ip;
IBmax=IB;
PaBar=a1a*Iamax+a2a*Iamax^2;
EB=0.95;
a1B=sigma*EB;
dnB=[0.25,0.15,2.5,0.5];
dnB=micro.*dnB;
a2B=sigma^2*sum(rhoN.*dnB)/aB;
PB=a1B*IB+a2B*IB^2;
Ep=0.95;
a1p=sigma*Ep;
dnp=[0.25,0.15,2.5,0.5];
dnp=micro.*dnp;
a2p=sigma^2*sum(rhoN.*dnp)/aP;
Pp=a1p*Ip+a2p*Ip^2;
Ppst=a1p*Ipst+a2p*Ipst^2;
PBst=PB;
P=2*PaBar+PB+PBst+Pp+Ppst;
d1=DeltaJB*k1*(1-tauB*(gammaN+g*Sssini));
d2=(gammaN+g*Sssini)*(Naini+k1*tauB*DeltaJB)-g*Sssini*Na0;
d3=a1a*d1+2*a2a*d1*d2;
d4=a2a*d1^2;
d5=a1a*d2+a2a*d2;
y01=(R2+R3)*P+T0;
y02=R3*P+T0;
y03=R1*P;
%%
m11=-1/(R2*Cc);
m12=-m11;
m13=1/(Cc*R1);
% m14=0;
m21=1/(R2*Cs);
m22=-(1/R2+1/R3)*(1/Cs);
m23=0;
% m24=0;
m31=6/(R2*Cc);
m32=-6/(R2*Cc);
m33=-6/(Cc*R1)-1/(R1*C);
% U1=1/Cc*(x*(d1*exp(-t/tauB)+d2-a1a/(2*a2a)+1/(2*a2a)*sqrt(a1a^2-4*a2a*d5+...
% 8*a2a*PaBar-4*a2a*d3*exp(-t/tauB)-4*a2a*d4*exp(-2*t/tauB))+Ires)^2)%Pc/Cc;
% U3=-6*1/Cc*(x*(d1*exp(-t/tauB)+d2-a1a/(2*a2a)+1/(2*a2a)*sqrt(a1a^2-4*a2a*d5+...
% 8*a2a*PaBar-4*a2a*d3*exp(-t/tauB)-4*a2a*d4*exp(-2*t/tauB))+Ires)^2);%
%% solve the system T0
tspan=[1e-11 1e-3];
y0=[y01 y02 y03 1e-11];
[t,y] = ode15s(@odefcn,tspan,y0);
%% ODE function
function dydt = odefcn(~,y)
dydt=zeros(4,1);
dydt(1)= m11*y(1)+ m12*y(2)+ m13*y(3)+1/Cc*(x*(d1*exp(-y(4)/tauB)+d2-a1a/(2*a2a)+1/(2*a2a)*sqrt(a1a^2-4*a2a*d5+...
8*a2a*PaBar-4*a2a*d3*exp(-y(4)/tauB)-4*a2a*d4*exp(-2*y(4)/tauB))+Ires)^2);
dydt(2)= m21*y(1)+ m22*y(2)+ m23*y(3);
dydt(3)= m31*y(1)+ m32*y(2)+ m33*y(3)+-6*1/Cc*(x*(d1*exp(-y(4)/tauB)+d2-a1a/(2*a2a)+1/(2*a2a)*sqrt(a1a^2-4*a2a*d5+...
8*a2a*PaBar-4*a2a*d3*exp(-y(4)/tauB)-4*a2a*d4*exp(-2*y(4)/tauB))+Ires)^2);
end
end

Sign in to comment.

Answers (0)

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by