MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: NDSolve and Plot

  • To: mathgroup at smc.vnet.net
  • Subject: [mg117931] Re: NDSolve and Plot
  • From: Dominic Woerner <dominic.woerner at mpi-hd.mpg.de>
  • Date: Tue, 5 Apr 2011 06:41:15 -0400 (EDT)

Clear[x, y, =CE=93, =CE=B8]
H = {{215, -104.1, 5.1, -4.3, 4.7, -15.1, -7.8},
          {-104.1, 220, 32.6, 7.1, 5.4, 8.3, 0.8},
          {5.1, 32.6, 0, -46.8, 1, -8.1, 5.1},
          {-4.3, 7.1, -46.8, 125, -70.7, -14.7, -61.5},
          {4.7, 5.4, 1, -70.7, 450, 89.7, -2.5},
          {-15.1, 8.3, -8.1, -14.7, 89.7, 330, 32.7},
          {-7.8, 0.8, 5.1, -61.5, -2.5, 32.7, 280}};

n = Length[H];
density = Table[rho[i, j][t], {i, n}, {j, n}];
rhoVec = Flatten[density];
rhoVecBare = Flatten[Table[rho[i, j], {i, n}, {j, n}]];
unitaryTerm = -I*(H.density - density.H);
trappingTerm =
  Table[-=CE=93*rho[i, j][t]*(KroneckerDelta[3, i] + KroneckerDelta[3, j]),{i, 1,
    n}, {j, 1, n}];

eomRHS = Flatten[unitaryTerm + trappingTerm];
eom = Table[D[rhoVec[[i]], t] == eomRHS[[i]], {i, 1, Length[rhoVec]}];
initialCondition =
  Table[rho[i, j][0] ==
     Cos[x]^2*KroneckerDelta[i, 1]*KroneckerDelta[j, 1] +
      Sin[x]^2*KroneckerDelta[i, 6]*KroneckerDelta[j, 6] +
      Sin[x]*Cos[x]*Exp[I*y]*KroneckerDelta[i, 1]*KroneckerDelta[j, 6] +
      Sin[x]*Cos[x]*Exp[-I*y]*KroneckerDelta[i, 6]*KroneckerDelta[j, 1], {i,
     1, n}, {j, 1, n}] // Flatten;

para = {=CE=93 -> 6.23, x -> =CF=80/4, y -> =CF=80*j/4, =CE=B8 -> 0};
tend = 100;
=CE=B7 = 1 - Sum[rho[i, i][t], {i, 1, n}];
(*sol=Partition[Flatten[Table[{=CF=80*k/4,=CF=80*m/100,Chop[=CE=B7/.NDSolve[{eom/.para/.j->k \
//.l->m,initialCondition/.para /.j->k \
//.l->m},rhoVecBare,{t,0,tend}]]},{k,-4,4},{m,-100,100}]],3];*)
sol = Table[
   Chop[=CE=B7 /.
     NDSolve[{eom /. para /. j -> k //. l -> m,
       initialCondition /. para /. j -> k //. l -> m},
      rhoVecBare, {t, 0, tend}, MaxSteps -> 10^10]], {k, -4, 4}];
(*ListPlot3D[sol]*)
LogLinearPlot[sol, {t, 0.01, tend},
 PlotRange -> {0,
   1}(*,PlotLegend->{"-=CF=80","-3/4=CF=80","-=CF=80/2","=CF=80-/4","0","=CF=80/4","=CF=80/2","3/4=CF=80","=CF=80"}*)]


>>>> I use NDSolve to find a solution to a set of first order differential =
equations and then I want to plot the results.
>>>>
>>>> My problem is, even if I solve the ODE till t==100, the lines in the p=
lot stop at around t==10.
>>>> There is no error and increasing MaxSteps doesn't help.
>>>>
>>>> Cheers,
>>>> dom
>>>>
>>>>
>>>
>>> Dominic,
>>>
>>> that is hard to analyze; could you provide some code?
>>>
>>> Oliver
>>>
>>


  • Prev by Date: Re: Manipulate, how to slowdown animation
  • Next by Date: Re: PlotLegend and Show
  • Previous by thread: Re: NDSolve and Plot
  • Next by thread: Trigonometry