| > |
w[c]:=0.01:delta:=5:Omega:=5:w[0]:=1:gamma_0:=1:k[b]:=1:T:=100:
nu := 2*Pi*k[b]*T:
k1 := 2*k[b]*T*w[c]^2
*
Sum((w[c]*exp(-w[c]*tau)-nu*n*exp(nu*n*tau))/(w[c]^2-(nu*n)^2),n=-infinity..+infinity):
F:=2*gamma_0*w[c]^2*exp(-w[c]*tau)*sin(delta/Omega*(sin(Omega*t)-sin(Omega*(t-tau)))+w[0]*tau):
G:=cos(delta/Omega*(sin(Omega*t)-sin(Omega*(t-tau)))+w[0]*tau)*k1:
eq:=diff(rho(t),t)+Int(G*rho(tau),tau=0..t)=1/2*Int(G-F,tau=0..t):
|
| > |
# Is this approximation reasonable???
G := cos(delta/Omega*(sin(Omega*t)-sin(Omega*(t-tau)))+w[0]*tau)*2*k[b]*T*w[c]^2
*
(w[c]*exp(-w[c]*tau)/(w[c]^2) + 2*Sum((w[c]*exp(-w[c]*tau)-nu*n*exp(nu*n*tau))/(w[c]^2-(nu*n)^2),n=1..5));
G := cos(delta/Omega*(sin(Omega*t)-sin(Omega*(t-tau)))+w[0]*tau)*2*k[b]*T*w[c]^2
*
(w[c]*exp(-w[c]*tau)/(w[c]^2) + 2*add((w[c]*exp(-w[c]*tau)-nu*n*exp(-2*Pi*n*k[b]*T*tau))/(w[c]^2-(nu*n)^2),n=1..5));
|
|
|
(1) |
| > |
# an important observation:
# with practically no loss of precision the rightmost term in G can be replaced by exp(4.605242862-0.1000007268e-1*tau)
op(-1, G);
z1 := eval(log(op(-1, G)), tau=0);
z2 := evalf(eval(log(op(-1, G)), tau=1000));
model := exp(z1+(z2-z1)/1000*tau);
# note that this model is consistent with the approximaton above for it is
# almost equal to the first term in op(-&, g)
expand(model), op(1, op(-1, G));
plots:-display(
plot(log(op(-1, G)), tau=0..1000, color=blue, thickness=5, transparency=0.5, legend=typeset('op(-1, G)')),
plot(log(model), tau=0..1000, color=red, legend=typeset('model')),
gridlines=true
);
plots:-display(
plot(model/op(-1, G), tau=0..1000, adaptive=true, color=red, legend=typeset('model'/'op(-1, G)')),
gridlines=true
);
|
| > |
# for saving computational time one could rewrite k1 this way:
'k1' = model;
|
|
|
(2) |
| > |
# A few more observations
# let's use "model" onstead of "op(-1, G)" in the expression of G,
# next define IGmodel(t) as the integral of G-F from 0 to t
G__model := t -> eval(cos(delta/Omega*(sin(Omega*t)-sin(Omega*(t-tau)))+w[0]*tau)*model)
|
|
![proc (t) options operator, arrow; eval(cos(delta*(sin(Omega*t)-sin(Omega*(t-tau)))/Omega+w[0]*tau)*model) end proc](/view.aspx?sf=274937_Answer/a9f38c21fa4e8888fd506e379b321e23.gif)
|
(3) |
| > |
F := unapply(F, t);
IG__model := t -> Int(G__model(t)-F(t), tau=0..t);
IG__model(t)
|
|
|
(4) |
| > |
# conjecture 1: limit(IG__model(t), t=+infinity)=0
# conjecture 2: the integral is periodic with perdiod 2*Pi
plot(evalf(eval(IG__model('t'))), t=0..100, gridlines=true, size=[1200, 400])
|
| > |
# Now we are going to build a numerical scheme to iteratively construct rho(t)
# To assess rho(t) at some time t=t[n+1] we start from its value rho(t[n]) at time t[n].
# When we have computed rho(t[n]) we knew the (at least theotitically) the value of IG__model(t[n]).
# The transition rho(t[n]) --> rho(t[n+1]) does not require computing IG_model(t[n+1]), but just
# the increment Int(G__model(t)-F, tau=t[n]..t[n+1]);
# This observation will save a lot of computational time
#
# thus let's define this quantity:
IG := (t, a, b) -> Int(G__model(t) , tau=a..b);
IGF := (t, a, b) -> Int(G__model(t)-F(t), tau=a..b);
|
| |
|
| |
| > |
Diff(rho(t),t)=-Int('G__model(t, tau)'*rho(tau),tau=0..t[n+1])+1/2*Int('G__model(t, tau)'-'F(t, tau)',tau=0..t[n+1]);
printf("\nAt time t=t[n+1] Diff(rho(t),t) can be approximated by\n");
Diff(rho(t),t)=(rho(t[n+1])-rho(t[n]))/dt;
printf("\nAt time t=t[n+1] the integral equation can be written this way\n");
rho(t[n+1])=rho(t[n])+dt*( -Int('G__model(t, tau)'*rho(tau),tau=0..t[n])-Int('G__model(t, tau)'*rho(tau),tau=t[n]..t[n+1])+1/2*Int('G__model(t, tau)'-'F(t, tau)',tau=0..t[n])+1/2*Int('G__model(t, tau)'-'F(t, tau)',tau=t[n]..t[n+1]) );
printf("\n...or, while introducing IGF and considering that the time t in IGF(t, t[n], t[n+1])\n");
printf("can be taken as (t[n]+t[n+1])/2 = t[n+1/2]\n");
rho(t[n+1])=rho(t[n])+dt*( -Int('G__model(t, tau)'*rho(tau),tau=0..t[n])-Int('G__model(t, tau)'*rho(tau),tau=t[n]..t[n+1])+1/2*'IGF(t, 0, t[n])'+1/2*'IGF(t[n+1/2], t[n], t[n+1])');
gh := `#mrow(mover(mrow(mo(G),mo("."),mo("ρ")),mo("^")))`:
printf("\nThis relation can be rewritten in this more concise form\n");
rho(t[n+1])=rho(t[n])+dt*( -gh(0..t[n])-gh(t[n]..t[n+1])+1/2*'IGF(t, 0, t[n])'+1/2*'IGF(t[n+1/2], t[n], t[n+1])') ;
printf("\nLet's remark that the two quantities below are known at time t=t[n+1] \n");
-gh(0..t[n]) and 1/2*'IGF(t, 0, t[n])';
printf("\nwhich means that the only ones to assess are\n");
-gh(t[n]..t[n+1]) and 1/2*'IGF(t, t[n], t[n+1])';
printf("\nThe second one poses no problem and can be computed by using evalf(eval('IGF(t[n+1/2], t[n], t[n+1])'))\n\n");
printf("and thus the only difficulty lies in the estimation of \n");
-gh(t[n]..t[n+1]);
printf("\n\nA simple way to do that is to consider that rho is constant over the interval t[n]..t[n+1]");
printf("A classical notation is write rho(t) this way\n");
rho(t) = rho(t[n+alpha]);
printf("where alpha is a number between 0 and 1 included:\n");
printf(" alpha = 0 leads to an explicit integration scheme\n");
printf(" alpha = 1 leads to an implicit integration scheme\n");
printf(" alpha = 1/2 will give the Crank-Nicholson scheme ot a multistep integration method\n");
printf("other choices are possible\n");
printf("\nExplicit schemes are usually only conditionnally stable, meaning the time step dt must not\n");
printf("exceed a maximum value to avoid spurious oscillations or numerical explosion\n");
printf("on the other way implicit schemes, if correctly designed, are unconditionnally stable\n");
printf("\nI use here the implicit (alpha=1) scheme\n");
rho(t[n+1])=rho(t[n])+dt*( -gh(0..t[n])-rho(t[n+1])*Int('G__model(t, tau)',tau=t[n]..t[n+1])+1/2*'IGF(t, 0, t[n])'+1/2*'IGF(t, t[n], t[n+1])' );
printf("\n... which can be rewritten as:\n");
rho(t[n+1])=rho(t[n])+dt*( -gh(0..t[n])-rho(t[n+1])*'IG(t[n], t[n+1])'+1/2*'IGF(t, 0, t[n])'+1/2*'IGF(t[n+1/2], t[n], t[n+1])' );
printf("\nAfter regrouping rho(t[n+1]) terms one gets\n");
rho(t[n+1])*(1+dt*'IG(t[n+1/2], t[n], t[n+1])')=rho(t[n])+dt*( -gh(0..t[n])+1/2*'IGF(t, 0, t[n])'+1/2*'IGF(t[n+1/2], t[n], t[n+1])' );
printf("\nAnd thus\n");
rho(t[n+1])=(rho(t[n])+dt*( -gh(0..t[n])+1/2*'IGF(t, 0, t[n])'+1/2*'IGF(t[n+1/2], t[n], t[n+1])' ))/(1+dt*'IG(t[n+1/2], t[n], t[n+1])');
printf("\n\nIn this expression the two terms\n");
'IGF(t[n+1/2], t[n], t[n+1])' and 'IG(t[n+1/2], t[n], t[n+1])';
printf("are computed using evalf(eval(""the term""))\n");
|
|
![Diff(rho(t), t) = -(Int(G__model(t, tau)*rho(tau), tau = 0 .. t[n+1]))+(1/2)*(Int(G__model(t, tau)-F(t, tau), tau = 0 .. t[n+1]))](/view.aspx?sf=274937_Answer/99a3e02d7e99dffd31c9b0405b1a6746.gif)
|
|
|
At time t=t[n+1] Diff(rho(t),t) can be approximated by
|
|
|
At time t=t[n+1] the integral equation can be written this way
|
|
|
...or, while introducing IGF and considering that the time t in IGF(t, t[n], t[n+1])
can be taken as (t[n]+t[n+1])/2 = t[n+1/2]
|
|
|
This relation can be rewritten in this more concise form
|
|
|
Let's remark that the two quantities below are known at time t=t[n+1]
|
|
|
which means that the only ones to assess are
|
|
|
The second one poses no problem and can be computed by using evalf(eval('IGF(t[n+1/2], t[n], t[n+1])'))
and thus the only difficulty lies in the estimation of
|
|
|
A simple way to do that is to consider that rho is constant over the interval t[n]..t[n+1]A classical notation is write rho(t) this way
|
|
|
where alpha is a number between 0 and 1 included:
alpha = 0 leads to an explicit integration scheme
alpha = 1 leads to an implicit integration scheme
alpha = 1/2 will give the Crank-Nicholson scheme ot a multistep integration method
other choices are possible
Explicit schemes are usually only conditionnally stable, meaning the time step dt must not
exceed a maximum value to avoid spurious oscillations or numerical explosion
on the other way implicit schemes, if correctly designed, are unconditionnally stable
I use here the implicit (alpha=1) scheme
|
|
|
... which can be rewritten as:
|
|
|
After regrouping rho(t[n+1]) terms one gets
|
|
|
In this expression the two terms
|
|
|
are computed using evalf(eval("the term"))
|
|
|
| > |
Dt := 0.01:
tend := 100:
times := [seq(0..tend, Dt)]:
nstep := numelems(times):
sol := Vector(nstep):
cumulative_Ghat := 0;
cumulative_GFhat := 0;
increment_Ghat := Dt * evalf(eval(IG (0, 0, Dt)));
increment_GFhat := (1/2) * Dt * evalf(eval(IGF(0, 0, Dt)));
sol[1] := 0.5:
|
|
|
(5) |
| > |
for m from 2 to nstep do
sol[m] := (sol[m-1] - Dt*cumulative_Ghat + Dt*cumulative_GFhat + Dt*increment_GFhat) / (1 + Dt*increment_Ghat);
cumulative_Ghat := cumulative_Ghat + increment_Ghat * (sol[m]+sol[m-1])/2;
cumulative_GFhat := cumulative_GFhat + increment_GFhat;
increment_Ghat := Dt * evalf(eval(IG ((m-1/2)*Dt, (m-1)*Dt, m*Dt)));
increment_GFhat := (1/2) * Dt * evalf(eval(IGF((m-1/2)*Dt, (m-1)*Dt, m*Dt)));
end do:
sol;
plot([seq([times[n], sol[n]], n=1..nstep)])
|
|
