Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- hh[x_, t_] := h[w[x[t], t]];
- ww[x_, t_] := w[x[t], t];
- eq0 = ((w == -2 Sqrt[1 + bb] + 2 Sqrt[h + bb] -
- Sqrt[bb] Log[(-Sqrt[bb] + Sqrt[1 + bb])/(
- Sqrt[bb] + Sqrt[1 + bb])] +
- Sqrt[bb] Log[(-Sqrt[bb] + Sqrt[h + bb])/(
- Sqrt[bb] + Sqrt[h + bb])]) /. {h -> hh[x, t],
- w -> ww[x, t]}) /. {h[w[x[t], t]] -> h[x, t],
- w[x[t], t] -> w[x, t]};
- eq1a = ((x[t] == (w + Sqrt[h + bb]) t) /. {h -> hh[x, t],
- w -> ww[x, t]}) /. {h[w[x[t], t]] -> h[x, t],
- w[x[t], t] -> w[x, t]};
- eq1b = ((x'[t] == w + Sqrt[h + bb]) /. {h -> hh[x, t],
- w -> ww[x, t]}) /. {h[w[x[t], t]] -> h[x, t],
- w[x[t], t] -> w[x, t]};
- eq2 = (ReleaseHold[
- Hold[(0 == D[w, t] + (w + Sqrt[h + bb]) D[w, x])] /. {h ->
- h[x, t], w -> w[x, t]}]);
- eq3 = (ReleaseHold[
- Hold[(0 == D[h, t] + w D[h, x] + h D[w, x])] /. {h -> h[x, t],
- w -> w[x, t]}]);
- eq4a = 0 == Limit[w[x, t], x -> Infinity];
- eq4b = 0 == Limit[w[x, t], x -> -Infinity];
- eq4c = 0 == Limit[w[x, t], Abs[x] -> -Infinity];
- eq5a = 1 == Limit[h[x, t], x -> Infinity];
- eq5b = 1 == Limit[h[x, t], x -> Infinity];
- eq5c = (ReleaseHold[
- Hold[(h == 1 + 1/2 Sech[x/20])] /. {h -> hh[x, 0],
- w -> ww[x, t]}] /. {x[0] -> x}) /. {h[w[x, 0]] -> h[x, 0]};
- eq5d = ReleaseHold[Hold[(Derivative[1, 0][h][0, t] == 0)]];
- xmax = 150.;
- tmax = 25.;
- soln = NDSolve[
- ({eq0, eq2, h[Infinity, t] == 1, h[-Infinity, t] == 1,
- w[Infinity, t] == 0, w[-Infinity, t] == 0, eq5c,
- eq5d} /. {bb -> 2.}),
- {w[x, t], h[x, t]},
- {t, 0, tmax},
- {x, -xmax, xmax},
- MaxStepSize -> 0.01, MaxSteps -> 10^6, SolveDelayed -> True
- ]
- (*
- eq5a=((ReleaseHold[Hold[(h==1)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->
- Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[
- w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
- eq5b=((ReleaseHold[Hold[(h==1)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->-
- Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[
- w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
- eq4a=((ReleaseHold[Hold[(w==0)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->
- Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[
- w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
- eq4b=((ReleaseHold[Hold[(w==0)]/.{h->hh[x,t],w->ww[x,t]}]/.{x[t]->-
- Infinity})/.{h[w[Infinity,t]]->h[0,t],h[w[-Infinity,t]]->h[0,t]})/.{h[
- w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
- eq2=(ReleaseHold[Hold[(0==D[w,t]+(w+Sqrt[h+bb])D[w,x])]/.{h->hh[x,t],
- w->ww[x,t]}])/.{h[w[x[t],t]]->h[x,t],w[x[t],t]->w[x,t]};
- eq3=(ReleaseHold[Hold[(0==D[h,t]+w D[h,x]+h
- D[w,x])]/.{h->hh[x,t],w->ww[x,t]}])/.{h[w[x[t],t]]->h[x,t],w[x[t],t]->
- w[x,t]};
- eq4a=ReleaseHold[Hold[(Limit[w,x->Infinity]==0)]/.{x->x[t],h->hh[x,t],
- w->ww[x,t]}];
- eq4b=ReleaseHold[Hold[(Limit[w,x->-Infinity]==0)]/.{x->x[t],h->hh[x,t]
- ,w->ww[x,t]}];
- eq5a=ReleaseHold[Hold[(Limit[h,x->Infinity]==1)]/.{x->x[t],h->hh[x,t],
- w->ww[x,t]}];
- eq5b=ReleaseHold[Hold[(Limit[h,x->-Infinity]==1)]/.{x->x[t],h->hh[x,t]
- ,w->ww[x,t]}];
- *)
- (*
- soln=NDSolve[
- ({eq0,h[Infinity,t]==1,h[-Infinity,t]==1,w[Infinity,t]==0,w[-Infinity,
- t]==0,eq5c,eq5d}/.{bb->2.}),
- {w[x,t],h[x,t]},
- {t,0,tmax},
- {x,-xmax,xmax},
- MaxStepSize->0.01,MaxSteps->10^6,SolveDelayed->True
- ];
- soln=NDSolve[
- ({eq0,eq4a,eq5a,eq5c,eq5d}/.{bb->2.}),
- {w[x,t],h[x,t]},
- {t,0,tmax},
- {x,-xmax,xmax},
- MaxStepSize->0.01,MaxSteps->10^6,SolveDelayed->True
- ];
- soln=NDSolve[
- {eq0,eq1a,eq1b,eq2,eq3,eq4c,eq5a,eq5c,eq5d},
- {w[x,t],h[x,t]},
- {t,0,tmax},
- {x,-xmax,xmax},
- MaxStepSize->0.01,MaxSteps->10^6,SolveDelayed->True
- ];
- soln = NDSolve[
- {eq0, eq1a, eq1b, eq2, Limit[w[x[t], t], x -> Infinity] == 0,
- Limit[w[x[t], t], x -> -Infinity] == 0,
- Limit[h[w[x[t], t]], x -> Infinity] == 1,
- Limit[h[w[x[t], t]], x -> -Infinity] == 1,
- Limit[h[w[x[t], t]], t -> 0] == 1 + 1/2 Sech[x[t]/20], x[0] == 0},
- {w[x[t], t], h[w[x[t], t]]},
- {t, 0, tmax},
- MaxStepSize -> 0.01, MaxSteps -> 10^6, SolveDelayed -> True
- ]
- *)
- tmax = 100.; (* Max. time step *)
- xmax = 500.; (* Max. spatial point *)
- dtmaxf = 0.0001; (* Max. fractional step *)
- beta = 2.; (* value for plasma beta *)
- xfac = 20; (* normalization factor for Sech *)
- xbounds = {-xmax/2, xmax};
- tbounds = {-1, tmax};
- (* Define initial condition functions *)
- hatt0[x_, xo_] := 1 + 1/2 Sech[x/xo];
- wofhb[hh_, bb_] := Sqrt[bb] + Sqrt[1 + bb] (-2 Sqrt[1 + bb] + 2 Sqrt[hh + bb] -
- Sqrt[bb] Log[(-Sqrt[bb] + Sqrt[1 + bb])/()] + Sqrt[bb]
- Log[(-Sqrt[bb] + Sqrt[hh + bb])/(Sqrt[bb] + Sqrt[hh + bb])]);
- (* Define partial differential equations *)
- 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;
- eq1 = D[h[x, t], t] + w[x, t] D[h[x, t], x] + h[x, t] D[w[x, t], x] == 0;
- eq2 = h[x, 0] == hatt0[x, xfac];
- eq3 = w[x, 0] == (wofhb[h, be]) /. {h -> (hatt0[x, xfac])};
- eq4a = h[xmax, t] == 1;
- eq4b = w[xmax, t] == 0;
- eqall = {eq0, eq1, eq2, eq3, eq4a, eq4b};
- (* Find solutions to partial differential equations as function of [x,t] *)
- soln = Block[{c = 0}, NDSolve[(eqall /. {be -> beta}),
- {w[x, t], h[x, t]},
- {x, xbounds[[1]], xbounds[[2]]},
- {t, tbounds[[1]], tbounds[[2]]},
- MaxSteps -> 10^6, SolveDelayed -> True,
- Method -> {"MethodOfLines", "TemporalVariable" -> t},
- MaxStepFraction -> dtmaxf, SolveDelayed -> True,
- StepMonitor :> If[Mod[c, 1000] == 0, Print["Step # = ", c]; c++, c++]]]
- (* Define time steps and colors for plots *)
- tsolns = {0, 10, 20, 30, 40, 50};
- cols = {Black, Purple, Blue, Green, Orange, Red};
- (* Define plot outputs *)
- plots = Table[Plot[((Evaluate[(h[x, t] /. Flatten[soln][[2]])]) /. {t -> tsolns[[i]]}),
- {x, -100, 180}, PlotRange -> All,
- PerformanceGoal -> "Quality",
- PlotStyle -> {AbsoluteThickness[1], cols[[i]]}],
- {i, 1, Length[tsolns], 1}];
- (* Show plots *)
- Show[plots, ImageSize -> {600, 600}, ImageMargins -> 0, AspectRatio -> Full]
- eqall = N[{eq0, eq1, eq2, eq3, eq4a, eq4b}];
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement