Problem with solving the inhomogeneous systems of ODE
4 views (last 30 days)
Show older comments
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
Walter Roberson
on 1 Aug 2020
Your t is a vector of values, which forces a number of other variables to be a vector of values, including U1 and U3. So in
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;
the first and the third of those wants to be a vector of values, one for each different t; meanwhile the second of the lines only wants to be a scalar.
As your U1 and U3 are built from "time", perhaps you should be using the original t vector together with the value passed in as the first parameter to odefcn to interpolate a particular U1 and U3 value ?
Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!