Prev: parse expression into expression tree
Next: find vector
From: gianluca messina on 19 May 2010 10:33
Hi,
i'm trying to solve this nonlinear differential equation system:

T' = (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*dirac(t-t0))./C;
Is' = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;

T and Is are function of time, like dirac(t-t0). I declare the system in a M-file where i give the value of the other parameters. Here ie the code:

function dy= sistema_non_omog (t,y)
T=y(1);
Is=y(2);
Rtes=Rn/2*(1+tanh(T-Tc+Is*1e-6*cost)./D);

dy = zeros(2,1);
dy(1) = (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*dirac(t-t0))./C;
dy(2) = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;
end

I recall this in my principle M-file with this code:

%solver
[t,y] = ode45(@sistema_non_omog,tspan,y0)

the file work, but when i plot T and Is, what i get is not what i expect. Maybe the problem is the dirac function? Any idea?

Thanks
From: Torsten Hennig on 19 May 2010 22:42
> Hi,
> i'm trying to solve this nonlinear differential
> equation system:
>
> T' =
> (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*dirac(t-t0))./C;
> Is' = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;
>
> T and Is are function of time, like dirac(t-t0). I
> declare the system in a M-file where i give the value
> of the other parameters. Here ie the code:
>
> function dy= sistema_non_omog (t,y)
> T=y(1);
> Is=y(2);
> Rtes=Rn/2*(1+tanh(T-Tc+Is*1e-6*cost)./D);
>
> dy = zeros(2,1);
> dy(1) =
> (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*dirac(t-t0))./C;
> dy(2) = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;
> end
>
> I recall this in my principle M-file with this code:
>
> %solver
> [t,y] = ode45(@sistema_non_omog,tspan,y0)
>
> the file work, but when i plot T and Is, what i get
> is not what i expect. Maybe the problem is the dirac
> function?

Yes, of course - how can you be sure the solver even
hits t0 exactly ? And what does a numerical method do
with a source term with value infinity ?

Choose a small intervall [t0-eps;t0+eps]
and distribute the mass 1 onto this intervall.

> Any idea?
>
> Thanks

Best wishes
Torsten.
From: gianluca messina on 20 May 2010 08:46
Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <1628349489.192154.1274337793968.JavaMail.root(a)gallium.mathforum.org>...
> > Hi,
> > i'm trying to solve this nonlinear differential
> > equation system:
> >
> > T' =
> > (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*dirac(t-t0))./C;
> > Is' = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;
> >
> > T and Is are function of time, like dirac(t-t0). I
> > declare the system in a M-file where i give the value
> > of the other parameters. Here ie the code:
> >
> > function dy= sistema_non_omog (t,y)
> > T=y(1);
> > Is=y(2);
> > Rtes=Rn/2*(1+tanh(T-Tc+Is*1e-6*cost)./D);
> >
> > dy = zeros(2,1);
> > dy(1) =
> > (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*dirac(t-t0))./C;
> > dy(2) = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;
> > end
> >
> > I recall this in my principle M-file with this code:
> >
> > %solver
> > [t,y] = ode45(@sistema_non_omog,tspan,y0)
> >
> > the file work, but when i plot T and Is, what i get
> > is not what i expect. Maybe the problem is the dirac
> > function?
>
> Yes, of course - how can you be sure the solver even
> hits t0 exactly ? And what does a numerical method do
> with a source term with value infinity ?
>
> Choose a small intervall [t0-eps;t0+eps]
> and distribute the mass 1 onto this intervall.
>
> > Any idea?
> >
> > Thanks
>
> Best wishes
> Torsten.

well,
the dirac is not what i want, because i'd like to perturb my system with a pulse narrow about 100 ns and high epsilon*E, so i've changed the system in:

T=y(1);
Is=y(2);

Rtes=Rn/2*(1+tanh(T-Tc+Is*1e-6*cost)./D);

dy = zeros(2,1);
dy(1) = (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*sinc(5000000*(t-t0)))./C;
dy(2) = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;

where i substitute dirac(t) with sinc(t).

But nothing is changed from before, moreover if i put E=0 (i.e. no pulse perturb the system) the result is the same.

I write the two file because i have no idea what the error is:

% file where i define the system
function dy= sistema_non_omog_1 (t,y)

global Rn Tc cost D

t0=5e-6;
epsilon=0.5;
E=0.0001516; % fJ, 1310 nm
k=18876086.4;
Rp = 7.5;% mOhm
Rsh = 9;% mOhm
Ts = 0.054; % K
Tc= 0.105; % K
Ib=9; % uA
C=15;% fJ/K
L=350;% nH

T=y(1);
Is=y(2);

% note that Rtes depend on T and IS
Rtes=Rn/2*(1+tanh(T-Tc+Is*1e-6*cost)./D);

dy = zeros(2,1);
dy(1) = (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*sinc(5000000*(t-t0)))./C;
dy(2) = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;
end
---------------------------------------------------------------------------------------------

% principle file
global t Rn Tc cost D k Rp Rsh

D = 0.0031111;
cost = 5101.8269;
Rn=1600; % mohm

% to be sure the solver hits the pulse
tspan = 0:5e-9:100e-6;

y0=[0.093173;1.793]; %[To(K);Is0(uA)]

%solver ode15s because i deal with a stiff problem (i guess for what i've understood about stiff problem)
[t,y] = ode15s(@sistema_non_omog_1,tspan,y0)
T1=y(:,1);
Is1=y(:,2);

Rtes1=Rn/2*(1+tanh((T1-Tc+Is1*1e-6*cost)./D)); %mOhm
Ptes1=Rtes1.*(Is1).^2; %fW

%graph
subplot(2,2,1)
plot (t,Is1,'.')
xlabel ('t (us)'),ylabel('I_s (uA)'),grid
ylim([0 max(Is1)*1.1])
title(['k=',num2str(k),' Rsh=',num2str(Rsh),' Rp=',num2str(Rp),' D=',num2str(D),' cost=',num2str(cost)])%'Ts=',num2str(Ts),'Tc=',num2str(Tc),
subplot(2,2,2)
plot (t,Rtes1,'.')
xlabel ('t (us)'),ylabel ('R_{tes} (ohm)'),grid
subplot(2,2,3)
plot (t,Ptes1,'.')
xlabel ('t (us)'),ylabel ('P_{tes} (fW)'),grid
subplot(2,2,4)
plot (t,T1,'.')
xlabel ('t (us)'),ylabel ('T (mK)'),grid
From: Torsten Hennig on 20 May 2010 05:23
> Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote
> in message
> <1628349489.192154.1274337793968.JavaMail.root(a)gallium
> .mathforum.org>...
> > > Hi,
> > > i'm trying to solve this nonlinear differential
> > > equation system:
> > >
> > > T' =
> > >
> (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*dirac(t-t0))./C;
> > > Is' = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;
> > >
> > > T and Is are function of time, like dirac(t-t0).
> I
> > > declare the system in a M-file where i give the
> value
> > > of the other parameters. Here ie the code:
> > >
> > > function dy= sistema_non_omog (t,y)
> > > T=y(1);
> > > Is=y(2);
> > > Rtes=Rn/2*(1+tanh(T-Tc+Is*1e-6*cost)./D);
> > >
> > > dy = zeros(2,1);
> > > dy(1) =
> > >
> (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*dirac(t-t0))./C;
> > > dy(2) = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;
> > > end
> > >
> > > I recall this in my principle M-file with this
> code:
> > >
> > > %solver
> > > [t,y] = ode45(@sistema_non_omog,tspan,y0)
> > >
> > > the file work, but when i plot T and Is, what i
> get
> > > is not what i expect. Maybe the problem is the
> dirac
> > > function?
> >
> > Yes, of course - how can you be sure the solver
> even
> > hits t0 exactly ? And what does a numerical method
> do
> > with a source term with value infinity ?
> >
> > Choose a small intervall [t0-eps;t0+eps]
> > and distribute the mass 1 onto this intervall.
> >
> > > Any idea?
> > >
> > > Thanks
> >
> > Best wishes
> > Torsten.
>
> well,
> the dirac is not what i want, because i'd like to
> perturb my system with a pulse narrow about 100 ns
> and high epsilon*E, so i've changed the system in:
>
> T=y(1);
> Is=y(2);
>
> Rtes=Rn/2*(1+tanh(T-Tc+Is*1e-6*cost)./D);
>
> dy = zeros(2,1);
> dy(1) =
> (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*sinc(5000000*(t-t
> 0)))./C;
> dy(2) = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;
>
> where i substitute dirac(t) with sinc(t).
>
> But nothing is changed from before, moreover if i put
> E=0 (i.e. no pulse perturb the system) the result is
> the same.
>
> I write the two file because i have no idea what the
> error is:
>
> % file where i define the system
> function dy= sistema_non_omog_1 (t,y)
>
> global Rn Tc cost D
>
> t0=5e-6;
> epsilon=0.5;
> E=0.0001516; % fJ, 1310 nm
> k=18876086.4;
> Rp = 7.5;% mOhm
> Rsh = 9;% mOhm
> Ts = 0.054; % K
> Tc= 0.105; % K
> Ib=9; % uA
> C=15;% fJ/K
> L=350;% nH
>
> T=y(1);
> Is=y(2);
>
> % note that Rtes depend on T and IS
> Rtes=Rn/2*(1+tanh(T-Tc+Is*1e-6*cost)./D);
>
> dy = zeros(2,1);
> dy(1) =
> (Rtes*Is.^2-k*(T.^5-Ts.^5)+epsilon*E*sinc(5000000*(t-t
> 0)))./C;
> dy(2) = (Ib*Rsh-Is.*(Rp+Rtes+Rsh))./L;
> end
> ------------------------------------------------------
> ---------------------------------------
>
> % principle file
> global t Rn Tc cost D k Rp Rsh
>
> D = 0.0031111;
> cost = 5101.8269;
> Rn=1600; % mohm
>
> % to be sure the solver hits the pulse
> tspan = 0:5e-9:100e-6;
>
> y0=[0.093173;1.793]; %[To(K);Is0(uA)]
>
> %solver ode15s because i deal with a stiff problem (i
> guess for what i've understood about stiff problem)
> [t,y] = ode15s(@sistema_non_omog_1,tspan,y0)
> T1=y(:,1);
> Is1=y(:,2);
>
> Rtes1=Rn/2*(1+tanh((T1-Tc+Is1*1e-6*cost)./D)); %mOhm
> Ptes1=Rtes1.*(Is1).^2; %fW
>
> %graph
> subplot(2,2,1)
> plot (t,Is1,'.')
> xlabel ('t (us)'),ylabel('I_s (uA)'),grid
> ylim([0 max(Is1)*1.1])
> title(['k=',num2str(k),' Rsh=',num2str(Rsh),'
> Rp=',num2str(Rp),' D=',num2str(D),'
> ,'
>
> cost=',num2str(cost)])%'Ts=',num2str(Ts),'Tc=',num2st
> r(Tc),
> subplot(2,2,2)
> plot (t,Rtes1,'.')
> xlabel ('t (us)'),ylabel ('R_{tes} (ohm)'),grid
> subplot(2,2,3)
> plot (t,Ptes1,'.')
> xlabel ('t (us)'),ylabel ('P_{tes} (fW)'),grid
> subplot(2,2,4)
> plot (t,T1,'.')
> xlabel ('t (us)'),ylabel ('T (mK)'),grid

Didn't you read my reply ?
The extra term only influences the solution at a
single time t=t0 which a numerical solver does
not recognize.
Distribute the extra term over a longer time interval.

Best wishes
Torsten.
 | 
Pages: 1
Prev: parse expression into expression tree
Next: find vector