Advertisement
Guest User

Untitled

a guest
Sep 22nd, 2017
66
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 5.80 KB | None | 0 0
  1. hh[x_, t_] := h[w[x[t], t]];
  2. ww[x_, t_] := w[x[t], t];
  3.  
  4. eq0 = ((w == -2 Sqrt[1 + bb] + 2 Sqrt[h + bb] -
  5. Sqrt[bb] Log[(-Sqrt[bb] + Sqrt[1 + bb])/(
  6. Sqrt[bb] + Sqrt[1 + bb])] +
  7. Sqrt[bb] Log[(-Sqrt[bb] + Sqrt[h + bb])/(
  8. Sqrt[bb] + Sqrt[h + bb])]) /. {h -> hh[x, t],
  9. w -> ww[x, t]}) /. {h[w[x[t], t]] -> h[x, t],
  10. w[x[t], t] -> w[x, t]};
  11. eq1a = ((x[t] == (w + Sqrt[h + bb]) t) /. {h -> hh[x, t],
  12. w -> ww[x, t]}) /. {h[w[x[t], t]] -> h[x, t],
  13. w[x[t], t] -> w[x, t]};
  14. eq1b = ((x'[t] == w + Sqrt[h + bb]) /. {h -> hh[x, t],
  15. w -> ww[x, t]}) /. {h[w[x[t], t]] -> h[x, t],
  16. w[x[t], t] -> w[x, t]};
  17. eq2 = (ReleaseHold[
  18. Hold[(0 == D[w, t] + (w + Sqrt[h + bb]) D[w, x])] /. {h ->
  19. h[x, t], w -> w[x, t]}]);
  20. eq3 = (ReleaseHold[
  21. Hold[(0 == D[h, t] + w D[h, x] + h D[w, x])] /. {h -> h[x, t],
  22. w -> w[x, t]}]);
  23. eq4a = 0 == Limit[w[x, t], x -> Infinity];
  24. eq4b = 0 == Limit[w[x, t], x -> -Infinity];
  25. eq4c = 0 == Limit[w[x, t], Abs[x] -> -Infinity];
  26. eq5a = 1 == Limit[h[x, t], x -> Infinity];
  27. eq5b = 1 == Limit[h[x, t], x -> Infinity];
  28. eq5c = (ReleaseHold[
  29. Hold[(h == 1 + 1/2 Sech[x/20])] /. {h -> hh[x, 0],
  30. w -> ww[x, t]}] /. {x[0] -> x}) /. {h[w[x, 0]] -> h[x, 0]};
  31. eq5d = ReleaseHold[Hold[(Derivative[1, 0][h][0, t] == 0)]];
  32.  
  33. xmax = 150.;
  34. tmax = 25.;
  35. soln = NDSolve[
  36. ({eq0, eq2, h[Infinity, t] == 1, h[-Infinity, t] == 1,
  37. w[Infinity, t] == 0, w[-Infinity, t] == 0, eq5c,
  38. eq5d} /. {bb -> 2.}),
  39. {w[x, t], h[x, t]},
  40. {t, 0, tmax},
  41. {x, -xmax, xmax},
  42. MaxStepSize -> 0.01, MaxSteps -> 10^6, SolveDelayed -> True
  43. ]
  44.  
  45. (*
  46. eq5a=((ReleaseHold[Hold[(h==1)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->
  47. Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[
  48. w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
  49. eq5b=((ReleaseHold[Hold[(h==1)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->-
  50. Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[
  51. w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
  52. eq4a=((ReleaseHold[Hold[(w==0)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->
  53. Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[
  54. w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
  55. eq4b=((ReleaseHold[Hold[(w==0)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->-
  56. Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[
  57. w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
  58. eq2=(ReleaseHold[Hold[(0==D[w,t]+(w+Sqrt[h+bb])D[w,x])]/.{h->hh[x,t],
  59. w->ww[x,t]}])/.{h[w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
  60. eq3=(ReleaseHold[Hold[(0==D[h,t]+w D[h,x]+h
  61. D[w,x])]/.{h->hh[x,t],w->ww[x,t]}])/.{h[w[x[t],t]]->h[x,t],w[x[t],t]->
  62. w[x,t]};
  63. eq4a=ReleaseHold[Hold[(Limit[w,x->Infinity]==0)]/.{x->x[t],h->hh[x,t],
  64. w->ww[x,t]}];
  65. eq4b=ReleaseHold[Hold[(Limit[w,x->-Infinity]==0)]/.{x->x[t],h->hh[x,t]
  66. ,w->ww[x,t]}];
  67. eq5a=ReleaseHold[Hold[(Limit[h,x->Infinity]==1)]/.{x->x[t],h->hh[x,t],
  68. w->ww[x,t]}];
  69. eq5b=ReleaseHold[Hold[(Limit[h,x->-Infinity]==1)]/.{x->x[t],h->hh[x,t]
  70. ,w->ww[x,t]}];
  71. *)
  72.  
  73. (*
  74. soln=NDSolve[
  75. ({eq0,h[Infinity,t]==1,h[-Infinity,t]==1,w[Infinity,t]==0,w[-Infinity,
  76. t]==0,eq5c,eq5d}/.{bb->2.}),
  77. {w[x,t],h[x,t]},
  78. {t,0,tmax},
  79. {x,-xmax,xmax},
  80. MaxStepSize->0.01,MaxSteps->10^6,SolveDelayed->True
  81. ];
  82. soln=NDSolve[
  83. ({eq0,eq4a,eq5a,eq5c,eq5d}/.{bb->2.}),
  84. {w[x,t],h[x,t]},
  85. {t,0,tmax},
  86. {x,-xmax,xmax},
  87. MaxStepSize->0.01,MaxSteps->10^6,SolveDelayed->True
  88. ];
  89. soln=NDSolve[
  90. {eq0,eq1a,eq1b,eq2,eq3,eq4c,eq5a,eq5c,eq5d},
  91. {w[x,t],h[x,t]},
  92. {t,0,tmax},
  93. {x,-xmax,xmax},
  94. MaxStepSize->0.01,MaxSteps->10^6,SolveDelayed->True
  95. ];
  96.  
  97. soln = NDSolve[
  98. {eq0, eq1a, eq1b, eq2, Limit[w[x[t], t], x -> Infinity] == 0,
  99. Limit[w[x[t], t], x -> -Infinity] == 0,
  100. Limit[h[w[x[t], t]], x -> Infinity] == 1,
  101. Limit[h[w[x[t], t]], x -> -Infinity] == 1,
  102. Limit[h[w[x[t], t]], t -> 0] == 1 + 1/2 Sech[x[t]/20], x[0] == 0},
  103. {w[x[t], t], h[w[x[t], t]]},
  104. {t, 0, tmax},
  105. MaxStepSize -> 0.01, MaxSteps -> 10^6, SolveDelayed -> True
  106. ]
  107. *)
  108.  
  109. tmax = 100.; (* Max. time step *)
  110. xmax = 500.; (* Max. spatial point *)
  111. dtmaxf = 0.0001; (* Max. fractional step *)
  112. beta = 2.; (* value for plasma beta *)
  113. xfac = 20; (* normalization factor for Sech *)
  114. xbounds = {-xmax/2, xmax};
  115. tbounds = {-1, tmax};
  116. (* Define initial condition functions *)
  117. hatt0[x_, xo_] := 1 + 1/2 Sech[x/xo];
  118. wofhb[hh_, bb_] := Sqrt[bb] + Sqrt[1 + bb] (-2 Sqrt[1 + bb] + 2 Sqrt[hh + bb] -
  119. Sqrt[bb] Log[(-Sqrt[bb] + Sqrt[1 + bb])/()] + Sqrt[bb]
  120. Log[(-Sqrt[bb] + Sqrt[hh + bb])/(Sqrt[bb] + Sqrt[hh + bb])]);
  121. (* Define partial differential equations *)
  122. eq0 = D[w[x, t], t] + w[x, t] D[w[x, t], x] + D[(h[x, t] + be Log[h[x, t]]), x] == 0;
  123. eq1 = D[h[x, t], t] + w[x, t] D[h[x, t], x] + h[x, t] D[w[x, t], x] == 0;
  124. eq2 = h[x, 0] == hatt0[x, xfac];
  125. eq3 = w[x, 0] == (wofhb[h, be]) /. {h -> (hatt0[x, xfac])};
  126. eq4a = h[xmax, t] == 1;
  127. eq4b = w[xmax, t] == 0;
  128. eqall = {eq0, eq1, eq2, eq3, eq4a, eq4b};
  129.  
  130. (* Find solutions to partial differential equations as function of [x,t] *)
  131. soln = Block[{c = 0}, NDSolve[(eqall /. {be -> beta}),
  132. {w[x, t], h[x, t]},
  133. {x, xbounds[[1]], xbounds[[2]]},
  134. {t, tbounds[[1]], tbounds[[2]]},
  135. MaxSteps -> 10^6, SolveDelayed -> True,
  136. Method -> {"MethodOfLines", "TemporalVariable" -> t},
  137. MaxStepFraction -> dtmaxf, SolveDelayed -> True,
  138. StepMonitor :> If[Mod[c, 1000] == 0, Print["Step # = ", c]; c++, c++]]]
  139.  
  140. (* Define time steps and colors for plots *)
  141. tsolns = {0, 10, 20, 30, 40, 50};
  142. cols = {Black, Purple, Blue, Green, Orange, Red};
  143. (* Define plot outputs *)
  144. plots = Table[Plot[((Evaluate[(h[x, t] /. Flatten[soln][[2]])]) /. {t -> tsolns[[i]]}),
  145. {x, -100, 180}, PlotRange -> All,
  146. PerformanceGoal -> "Quality",
  147. PlotStyle -> {AbsoluteThickness[1], cols[[i]]}],
  148. {i, 1, Length[tsolns], 1}];
  149. (* Show plots *)
  150. Show[plots, ImageSize -> {600, 600}, ImageMargins -> 0, AspectRatio -> Full]
  151.  
  152. eqall = N[{eq0, eq1, eq2, eq3, eq4a, eq4b}];
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement