Mathematica notebook for manuscript
Klika V: Pattern formation revisited within non-equilibrium thermodynamics: Burgers’ type equation
submitted to Biological Cybernetics 2021

NN=100;

(*randIC = RandomReal[UniformDistribution[{-0.01, 0.01}], NN+1];*)randIC={0.008468322241858513`,-0.0038596420569447197`,0.008684018261921854`,0.0005485710357122538`,-0.0029081904393727168`,-0.0017823763335907856`,0.002244229185843806`,-0.008111168050964148`,-0.0030479365932791803`,-0.007619572947200268`,-0.008117294470189607`,-0.006574465015407791`,-0.002430081583695353`,-0.007731994473206671`,0.004909768803176534`,0.009644215675422271`,-0.0011226446055206847`,-0.001323932240144203`,0.0009759807480402942`,-0.00027869387289218214`,-0.007511405180100399`,0.007740726843651395`,-0.009913751652711612`,-0.00028882272473379836`,-0.0015317882579189042`,0.005087387792514442`,-0.0012543229932629324`,-0.007889681743632762`,0.005762648579780122`,0.002105240351657621`,-0.0070394172877857805`,-0.007305984200420352`,-0.0009106220923104473`,-0.009972598201472022`,-0.006420133724939506`,0.007876014595470271`,0.0051369173131669`,0.008543002895672197`,-0.0020192675333090965`,0.008149889225058486`,-0.003835684661223644`,0.009797409095060716`,0.005944316385770317`,0.006826001790665226`,0.005902397647162105`,0.0056022017921813425`,-0.006556066985842376`,-0.009450429078459166`,0.0021680537991390914`,0.00023905267471401384`,0.009118065257567795`,0.009756478443435856`,-0.0016741906419891726`,0.008310243217345574`,0.0006668281922656387`,-0.0002552182937378099`,-0.0057933117892314985`,-0.0027564369239762992`,-0.0002045819725730695`,-0.00940298369147349`,-0.006019688226263208`,0.006593232437555614`,0.006794167519888469`,0.001792231621492113`,-0.00032373805242035675`,0.009738371524471219`,0.009613679750365416`,-0.00458292117977761`,0.006903211986609203`,0.007195808198647756`,0.000025466374971722106`,0.00021064740089821524`,0.006875909176554437`,0.00522582711240268`,-0.006230218039499011`,-0.009632243804867189`,0.0040424720079230425`,-0.003343820876719604`,0.00629183129022029`,-0.003258112436123558`,0.006318377403628284`,-0.0032530733380565814`,-0.002882597716339575`,0.008991945666113825`,0.0029352680917725665`,-0.006761182166516371`,0.00858859904337056`,0.005778857513179789`,0.0038519150225514373`,-0.005948074141475596`,0.004620043461836566`,-0.0011314433503448185`,-0.0020103053568665125`,-0.0008172639736635581`,0.00691496282438428`,-0.004384591052711138`,0.0014947864209754955`,0.005121565817013507`,-0.009506601928939686`,0.008484252739160779`,-0.009714062619114582`};

Fig1,6 Burger’s eq with a linear damping; Dirichlet BCs - evolution to a seemingly stationary state

L=100;NN=100;
ν=1;r= 1 (10^0);

iniCondE =  ListInterpolation[   Flatten[Prepend[Flatten[Append[ randIC[[ 1 ;; NN + 1]], {randIC[[NN + 1]]}]], {randIC[[ 1]]}]], {{0, L}}, InterpolationOrder->1];
Plot[iniCondE[x],{x,0,L}]
Plot[D[iniCondE[y],y]/.y->x,{x,0,L}]

iniCond = {u[0,x]==D[iniCondE[x],x]};
BC = {u[t, 0] ==0, u[t, L] == 0};
DEqns = {D[u[t,x],t] == ν D[u[t,x],{x,2}]-2 (10^0) u[t,x] D[u[t,x],x]+r u[t,x]};

τ =2 (10^1);

solNBurgers = NDSolve[Flatten[{DEqns,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.01, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->20000}}}]

(*test*)

Plot[{(u[t,0]/.solNBurgers),(u[t,L]/.solNBurgers)},{t,0,τ},PlotRange->All]

Plot[Evaluate[u[t, x] /. solNBurgers/. t -> 10], {x, 0, L},PlotStyle->{Thick,Red},PlotRange->All,PlotPoints->100,AxesLabel->{Style["x",20],Style["y",20]},TicksStyle->Directive[FontSize->18]]

Plot[Evaluate[u[t, x] /. solNBurgers/. t -> 12], {x, 0, L},PlotStyle->{Thick,Red},PlotRange->All,PlotPoints->100]

Plot[Evaluate[u[t, x] /. solNBurgers/. t -> 13], {x, 0, L},PlotStyle->{Thick,Red},PlotRange->All,PlotPoints->100]

Plot[Evaluate[u[t, x] /. solNBurgers/. t -> τ], {x, 0, L},PlotStyle->{Thick,Red},PlotRange->All,PlotPoints->100]

Plot3D[Evaluate[u[t, x] /. solNBurgers], {x, 0, L}, {t, 0, τ},PlotPoints->200,Mesh->None,PlotRange->All]

Show[Plot3D[Evaluate[u[t, x] /. solNBurgers], {x, 0, L}, {t, 0, τ},PlotPoints->200,Mesh->None,PlotRange->All,AxesLabel->{Style["x",20],Style["t",20],Style["y",20]},TicksStyle->Directive[FontSize->18]],ContourPlot3D[t==10,{x,0,L},{t,0,τ},{y,0,50},ContourStyle->Opacity[0.2]]]

Export["Fig1a_alt.pdf",%]

DensityPlot[Evaluate[u[t, x] /. solNBurgers],  {x, 0, L},{t, 0, τ}(*,PlotRange->All*)(*,ColorFunction->Hue*)(*"SunsetColors"*),LabelStyle->Directive[Medium],(*FrameLabel->{t, x},*)PlotPoints->400,PlotRange->All,FrameLabel->{Style["x",20],Style["t",20]},FrameTicksStyle->Directive[FontSize->18]]

τ=20;

DEcheck=Table[Evaluate[DEqns[[1,1]]-DEqns[[1,2]]/.solNBurgers[[1]]/.t->τ],{x,0,L,L/(2 (10^4))}];

ListPlot[DEcheck,AxesLabel->{Style["Nodes",20],Style["Error",20]},TicksStyle->Directive[FontSize->18]]
ListPlot[DEcheck,PlotRange->All,AxesLabel->{Style["Nodes",20],Style["Error",20]},TicksStyle->Directive[FontSize->18]]

First[DEcheck]
Last[DEcheck]

Sort[DEcheck,Abs[#2]<Abs[#1]&]

timeDerivativeCheck=Table[Evaluate[DEqns[[1,1]]/.solNBurgers[[1]]/.t->τ],{x,0,L,L/(5 (10^4))}];

ListPlot[timeDerivativeCheck]
ListPlot[timeDerivativeCheck,PlotRange->All,AxesLabel->{Style["Nodes",20],Style["Error",20]},TicksStyle->Directive[FontSize->18]]
Sort[timeDerivativeCheck,Abs[#2]<Abs[#1]&]

Fig 2- Stationary Solution; various BCs

NDSolve[{u[0]==0,u[2]==0,0==u''[x]+(((u'[x]))^2)+u[x] },u[x],{x,0,2}, Method->"BoundaryValues"->{"Shooting", "StartingInitialConditions"->{u[0] ==0, u'[0] ==10}}(*,Method->"StiffnessSwitching",AccuracyGoal->10,PrecisionGoal->10*)]
Plot[u[x]/.%,{x,0,2},PlotRange->All,PlotStyle->Thick,PlotTheme->"Monochrome",AxesLabel->{Style["x",20],Style["ρ",20]},TicksStyle->Directive[FontSize->18]]

Export["Fig2.pdf",%]

NDSolve[{u[0]==1,u[2]==1,0==u''[x]+(((u'[x]))^2)+(u[x]+1)+(u[x]-1) (u[x]-2)(u[x]-3)},u[x],{x,0,2}, Method->"BoundaryValues"->{"Shooting", "StartingInitialConditions"->{u[0] == 1, u'[0] ==1}}(*,Method->"StiffnessSwitching",AccuracyGoal->10,PrecisionGoal->10*)]
Plot[u[x]/.%,{x,0,2},PlotRange->All,PlotStyle->Thick,PlotTheme->"Monochrome"]

NDSolve[{u[0]==0,u'[10]==0,0==u''[x]+(10^0) (((u'[x]))^2)+(10^0) (u[x]-1)},u[x],{x,0,10}]
Plot[u[x]/.%,{x,0,10},PlotRange->All,PlotStyle->Thick,PlotTheme->"Monochrome"]

NDSolve[{u[0]==0,u[10]==0,0==u''[x]+(10^0) (((u'[x]))^2)+(10^0) (u[x]-1)},u[x],{x,0,10}]
Plot[u[x]/.%,{x,0,10},PlotRange->All,PlotStyle->Thick,PlotTheme->"Monochrome"]

NDSolve[{u[0]==0,u'[100]==1,0==u''[x]+(10^0) (((u'[x]))^2)+(10^0) (u[x]-1)},u[x],{x,0,100}]
Plot[u[x]/.%,{x,0,100},PlotRange->All,PlotStyle->Thick,PlotTheme->"Monochrome"]

Fig 3 - cubic reaction with Zero BCs; stab->instab transition

L =10;ν=(10^(-1));r=(10^(-1));

iniCondE =  ListInterpolation[   Flatten[Prepend[Flatten[Append[ randIC[[ 1 ;; 100 + 1]], {randIC[[100 + 1]]}]], {randIC[[ 1]]}]], {{0, L}}, InterpolationOrder->1];(*BC = zero flux*)
iniCond = {u[0,x]==0.995+iniCondE[x]};
BC = {Derivative[0,1][u][t, 0] ==0, Derivative[0,1][u][t, L] == 0};
DEqns = {D[u[t,x],t] == ν D[u[t,x],{x,2}]+κ(((D[u[t,x],x]))^2)+ r (-4 u[t,x]+4 ((u[t,x])^3))};

τ =50;

Plot[iniCond[[1,2]],{x,0,L}]

datapoints={};
For[k=0,k≤33,k++,solNRDnKED = NDSolve[Flatten[{DEqns/.κ->k,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}, StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
AppendTo[datapoints,{k,NIntegrate[(1/L)(((u[t,x]))^2)/.solNRDnKED/.t->τ,{x,0,L},Method->{Automatic,"SymbolicProcessing"->False}][[1]]}];Print[k,",",datapoints];
]

τ=50;
solNRDnKED = NDSolve[Flatten[{DEqns/.κ->31,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}, StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->100000}}}]

Plot[Evaluate[u[t, x] /. solNRDnKED/. t -> τ], {x, 0, L},PlotStyle->{Thick,Red}]

Plot3D[Evaluate[u[t, x] /. solNRDnKED], {x, 0, L}, {t, 0, τ},PlotPoints->50,PlotRange->All]

τ=15;
solNRDnKED = NDSolve[Flatten[{DEqns/.κ->32,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}, StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}]

Plot[Evaluate[u[t, x] /. solNRDnKED/. t -> 14.8], {x, 0, L},PlotStyle->{Thick,Red}]

Plot3D[Evaluate[u[t, x] /. solNRDnKED], {x, 0, L}, {t, 0, 14.8},PlotPoints->50,PlotRange->All]

Plot3D[Log[Evaluate[u[t, x] /. solNRDnKED]], {x, 0, L}, {t, 0, 14.8},PlotPoints->50,PlotRange->All]

τ=11;
solNRDnKED = NDSolve[Flatten[{DEqns/.κ->32,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}, StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->150000}}}]

Plot[Evaluate[u[t, x] /. solNRDnKED/. t -> 10.8], {x, 0, L},PlotStyle->{Thick,Red}]

Plot3D[Evaluate[u[t, x] /. solNRDnKED], {x, 0, L}, {t, 0, 10.8},PlotPoints->50,PlotRange->All]

τ=7;

For[k=33,k≤40,k++,solNRDnKED = NDSolve[Flatten[{DEqns/.κ->k,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}, StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
AppendTo[datapoints,{k,NIntegrate[(1/L)(((u[t,x]))^2)/.solNRDnKED/.t->τ,{x,0,L},Method->{Automatic,"SymbolicProcessing"->False}][[1]]}];Print[k];
]

ListPlot[datapoints,PlotRange->All,AxesLabel->{Style["K",20],Style["∥u(Subscript[∥, (L^2)])",20]},TicksStyle->Directive[FontSize->18]]

Export["Fig3a.pdf",%]

increasing the IC

L =10;
randIC=Get["ridge2"];

ν=(10^(-1));r=(10^(-1));

iniCondE =  ListInterpolation[   Flatten[Prepend[Flatten[Append[ randIC[[ 1 ;; 100 + 1]], {randIC[[100 + 1]]}]], {randIC[[ 1]]}]], {{0, L}}, InterpolationOrder->1];(*BC = zero flux*)
iniCond = {u[0,x]==0.995+1.5 iniCondE[x]};
BC = {Derivative[0,1][u][t, 0] ==0, Derivative[0,1][u][t, L] == 0};
DEqns = {D[u[t,x],t] == ν D[u[t,x],{x,2}]+κ(((D[u[t,x],x]))^2)+ r (-4 u[t,x]+4 ((u[t,x])^3))};

τ =50;

Plot[iniCond[[1,2]],{x,0,L}]

datapoints={};
For[k=0,k≤9,k++,solNRDnKED = NDSolve[Flatten[{DEqns/.κ->k,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}, StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
AppendTo[datapoints,{k,NIntegrate[(1/L)(((u[t,x]))^2)/.solNRDnKED/.t->τ,{x,0,L},Method->{Automatic,"SymbolicProcessing"->False}][[1]]}];Print[k];
]

datapoints

datapoints=datapoints[[;;8]]

τ=7;

For[k=8,k≤20,k++,solNRDnKED = NDSolve[Flatten[{DEqns/.κ->k,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}, StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
AppendTo[datapoints,{k,NIntegrate[(1/L)(((u[t,x]))^2)/.solNRDnKED/.t->τ,{x,0,L},Method->{Automatic,"SymbolicProcessing"->False}][[1]]}];Print[k];
]

ListPlot[datapoints[[;;15]],PlotRange->{{0,40},Automatic},AxesLabel->{Style["K",20],Style["∥u(Subscript[∥, (L^2)])",20]},TicksStyle->Directive[FontSize->18]]

Export["Fig3b.pdf",%]

Fig 4,5 - The study of the significance of K/D

Fig 4 - r=1
lin kin, Burgers eqn

ν=1;

L=100;r=1;

iniCondE =  ListInterpolation[   Flatten[Prepend[Flatten[Append[ randIC[[ 1 ;; 100 + 1]], {randIC[[100 + 1]]}]], {randIC[[ 1]]}]], {{0, L}}, InterpolationOrder->1];

iniCond = {u[0,x]==0.5D[iniCondE[x],x]};
BC = {u[t, 0] ==0, u[t, L] == 0};
DEqns = {D[u[t,x],t] == ν D[u[t,x],{x,2}]-2 κ u[t,x] D[u[t,x],x]+r u[t,x]};
τ=10;

datapoints={};For[k=0,k≤2,k=k+1/20,solNBurgers =  NDSolve[Flatten[{DEqns/.κ->k,iniCond,BC}], {u}, {t,0,τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.01, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
AppendTo[datapoints,{k,Max[Table[Abs[D[u[τ,y],y,y]]/.solNBurgers/.y->x,{x,0,L,L/(10^4)}]](*,Max[Table[Abs[D[u[τ,y],y,y,y]]/.solNBurgers/.y->x,{x,0,L,L/(10^4)}]]*),NIntegrate[(1/L)(((u[t,x]))^2)/.solNBurgers/.t->τ,{x,0,L},Method->{Automatic,"SymbolicProcessing"->False}][[1]],NIntegrate[(1/L)(((D[u[t,x],x,x]))^2)/.solNBurgers/.t->τ,{x,0,L},Method->{Automatic,"SymbolicProcessing"->False}][[1]]}];
Print[k,datapoints[[20 k+1]]];]

ListPlot[Table[{datapoints[[k,1]],(1/(datapoints[[1,4]]/(((datapoints[[1,2]]))^2)))datapoints[[k,4]]/(((datapoints[[k,2]]))^2)},{k,1,Length[datapoints]}]]
ListPlot[Table[{datapoints[[k,1]],(1/(datapoints[[1,4]]/(((datapoints[[1,2]]))^2)))datapoints[[k,4]]/(((datapoints[[k,2]]))^2)},{k,1,Length[datapoints]}],PlotRange->All]
ListPlot[Take[datapoints,All,{1,2}]/datapoints[[1,2]]]

ListPlot[Take[datapoints,All,{1,4,3}]/datapoints[[1,4]]]

solNRDnKED1 = NDSolve[Flatten[{DEqns/.κ->0.3,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];

solNRDnKED1 = NDSolve[Flatten[{DEqns/.κ->0.2,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
solNRDnKED2= NDSolve[Flatten[{DEqns/.κ->1,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
solNRDnKED3= NDSolve[Flatten[{DEqns/.κ->2,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];

solNRD = NDSolve[Flatten[{DEqns/.κ->0,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];

Plot[{Evaluate[u[t, x] /. solNRD/. t -> τ],Evaluate[u[t, x] /. solNRDnKED1/. t -> τ],Evaluate[u[t, x] /. solNRDnKED2/. t -> τ],Evaluate[u[t, x] /. solNRDnKED3/. t -> τ]}, {x, 0, L}(*,PlotStyle->{{Thick,Red},Blue}*)]
Plot[{Evaluate[u[t, x] /. solNRD/. t -> τ],Evaluate[u[t, x] /. solNRDnKED1/. t -> τ],Evaluate[u[t, x] /. solNRDnKED2/. t -> τ],Evaluate[u[t, x] /. solNRDnKED3/. t -> τ]}, {x, 0, L},PlotStyle->{Thick,Thick,Thick,Thick},AxesLabel-> {Style["x",20],Style["y",20]},TicksStyle->Directive[FontSize->18],PlotTheme->"Monochrome"]

Export["Fig5b.pdf",%]

Fig 4 - r=0.2
lin kin, Burgers eqn

ν=1;

L=100;r=0.2;

iniCondE =  ListInterpolation[   Flatten[Prepend[Flatten[Append[ randIC[[ 1 ;; 100 + 1]], {randIC[[100 + 1]]}]], {randIC[[ 1]]}]], {{0, L}}, InterpolationOrder->1];

iniCond = {u[0,x]==0.5D[iniCondE[x],x]};
BC = {u[t, 0] ==0, u[t, L] == 0};
DEqns = {D[u[t,x],t] == ν D[u[t,x],{x,2}]-2 κ u[t,x] D[u[t,x],x]+r u[t,x]};
τ=50;

datapoints={};For[k=0,k≤2,k=k+1/20,solNBurgers =  NDSolve[Flatten[{DEqns/.κ->k,iniCond,BC}], {u}, {t,0,τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.01, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
AppendTo[datapoints,{k,Max[Table[Abs[D[u[τ,y],y,y]]/.solNBurgers/.y->x,{x,0,L,L/(10^4)}]](*,Max[Table[Abs[D[u[τ,y],y,y,y]]/.solNBurgers/.y->x,{x,0,L,L/(10^4)}]]*),NIntegrate[(1/L)(((u[t,x]))^2)/.solNBurgers/.t->τ,{x,0,L},Method->{Automatic,"SymbolicProcessing"->False}][[1]],NIntegrate[(1/L)(((D[u[t,x],x,x]))^2)/.solNBurgers/.t->τ,{x,0,L},Method->{Automatic,"SymbolicProcessing"->False}][[1]]}];
Print[k,datapoints[[20 k+1]]];]

ListPlot[Table[{datapoints[[k,1]],(1/(datapoints[[1,4]]/(((datapoints[[1,2]]))^2)))datapoints[[k,4]]/(((datapoints[[k,2]]))^2)},{k,1,Length[datapoints]}]]
ListPlot[Table[{datapoints[[k,1]],(1/(datapoints[[1,4]]/(((datapoints[[1,2]]))^2)))datapoints[[k,4]]/(((datapoints[[k,2]]))^2)},{k,1,Length[datapoints]}],PlotRange->All]
ListPlot[Table[{datapoints[[k,1]],datapoints[[k,2]]/datapoints[[1,2]]},{k,1,Length[datapoints]}]]

ListPlot[Table[{datapoints[[k,1]],datapoints[[k,4]]/datapoints[[1,4]]},{k,1,Length[datapoints]}]]

solNRDnKED1 = NDSolve[Flatten[{DEqns/.κ->0.2,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
solNRDnKED2= NDSolve[Flatten[{DEqns/.κ->1,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
solNRDnKED3= NDSolve[Flatten[{DEqns/.κ->2,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];

solNRD = NDSolve[Flatten[{DEqns/.κ->0,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];

Plot[{Evaluate[u[t, x] /. solNRD/. t -> τ],Evaluate[u[t, x] /. solNRDnKED1/. t -> τ],Evaluate[u[t, x] /. solNRDnKED2/. t -> τ],Evaluate[u[t, x] /. solNRDnKED3/. t -> τ]}, {x, 0, L}(*,PlotStyle->{{Thick,Red},Blue}*)]
Plot[{Evaluate[u[t, x] /. solNRD/. t -> τ],Evaluate[u[t, x] /. solNRDnKED1/. t -> τ],Evaluate[u[t, x] /. solNRDnKED2/. t -> τ],Evaluate[u[t, x] /. solNRDnKED3/. t -> τ]}, {x, 0, L},PlotStyle->{Thick,Thick,Thick,Thick},AxesLabel-> {Style["x",20],Style["y",20]},TicksStyle->Directive[FontSize->18],PlotTheme->"Monochrome"]

Export["Fig5a.pdf",%]

datapoints02=datapoints;

Fig 4 - r=0.6
lin kin, Burgers eqn

ν=1;

L=100;r=0.6;

iniCondE =  ListInterpolation[   Flatten[Prepend[Flatten[Append[ randIC[[ 1 ;; 100 + 1]], {randIC[[100 + 1]]}]], {randIC[[ 1]]}]], {{0, L}}, InterpolationOrder->1];

iniCond = {u[0,x]==0.5D[iniCondE[x],x]};
BC = {u[t, 0] ==0, u[t, L] == 0};
DEqns = {D[u[t,x],t] == ν D[u[t,x],{x,2}]-2 κ u[t,x] D[u[t,x],x]+r u[t,x]};
τ=10/0.6;

datapoints={};For[k=0,k≤2,k=k+1/20,solNBurgers =  NDSolve[Flatten[{DEqns/.κ->k,iniCond,BC}], {u}, {t,0,τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.01, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
AppendTo[datapoints,{k,Max[Table[Abs[D[u[τ,y],y,y]]/.solNBurgers/.y->x,{x,0,L,L/(10^4)}]](*,Max[Table[Abs[D[u[τ,y],y,y,y]]/.solNBurgers/.y->x,{x,0,L,L/(10^4)}]]*),NIntegrate[(1/L)(((u[t,x]))^2)/.solNBurgers/.t->τ,{x,0,L},Method->{Automatic,"SymbolicProcessing"->False}][[1]],NIntegrate[(1/L)(((D[u[t,x],x,x]))^2)/.solNBurgers/.t->τ,{x,0,L},Method->{Automatic,"SymbolicProcessing"->False}][[1]]}];
Print[k,datapoints[[20 k+1]]];]

ListPlot[Table[{datapoints[[k,1]],(1/(datapoints[[1,4]]/(((datapoints[[1,2]]))^2)))datapoints[[k,4]]/(((datapoints[[k,2]]))^2)},{k,1,Length[datapoints]}]]
ListPlot[Table[{datapoints[[k,1]],(1/(datapoints[[1,4]]/(((datapoints[[1,2]]))^2)))datapoints[[k,4]]/(((datapoints[[k,2]]))^2)},{k,1,Length[datapoints]}],PlotRange->All]
ListPlot[Table[{datapoints[[k,1]],datapoints[[k,2]]/datapoints[[1,2]]},{k,1,Length[datapoints]}]]

ListPlot[Table[{datapoints[[k,1]],datapoints[[k,4]]/datapoints[[1,4]]},{k,1,Length[datapoints]}]]

solNRDnKED1 = NDSolve[Flatten[{DEqns/.κ->0.2,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
solNRDnKED2= NDSolve[Flatten[{DEqns/.κ->1,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];
solNRDnKED3= NDSolve[Flatten[{DEqns/.κ->2,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];

solNRD = NDSolve[Flatten[{DEqns/.κ->0,iniCond,BC}], {u}, {t,0, τ}, {x, 0, L}(*,MaxStepSize->0.1*)(*,AccuracyGoal->6,PrecisionGoal->7*)(*,MaxSteps->{100,(10^5)}*), StartingStepSize ->0.0001, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000,"MaxPoints"->50000}}}];

Plot[{Evaluate[u[t, x] /. solNRD/. t -> τ],Evaluate[u[t, x] /. solNRDnKED1/. t -> τ],Evaluate[u[t, x] /. solNRDnKED2/. t -> τ],Evaluate[u[t, x] /. solNRDnKED3/. t -> τ]}, {x, 0, L}(*,PlotStyle->{{Thick,Red},Blue}*)]
Plot[{Evaluate[u[t, x] /. solNRD/. t -> τ],Evaluate[u[t, x] /. solNRDnKED1/. t -> τ],Evaluate[u[t, x] /. solNRDnKED2/. t -> τ],Evaluate[u[t, x] /. solNRDnKED3/. t -> τ]}, {x, 0, L},PlotStyle->{Thick,Thick},PlotTheme->"Monochrome"]

datapoints06=datapoints;

together

ListPlot[{Table[{datapoints02[[k,1]],(1/(datapoints02[[1,4]]/(((datapoints02[[1,2]]))^2)))datapoints02[[k,4]]/(((datapoints02[[k,2]]))^2)},{k,1,Length[datapoints02]}],Table[{datapoints06[[k,1]],(1/(datapoints06[[1,4]]/(((datapoints06[[1,2]]))^2)))datapoints06[[k,4]]/(((datapoints06[[k,2]]))^2)},{k,1,Length[datapoints06]}],Table[{datapoints10[[k,1]],(1/(datapoints10[[1,4]]/(((datapoints10[[1,2]]))^2)))datapoints10[[k,4]]/(((datapoints10[[k,2]]))^2)},{k,1,Length[datapoints10]}]},PlotRange->All,PlotStyle->{Thick,Thick,Thick},PlotTheme->"Monochrome"(*PlotLegends->{0.2,0.6,1.0}*),AxesLabel->{Style["K",20],Style["(Subscript[κ, 2])",20]},TicksStyle->Directive[FontSize->18]]
Export["Fig4b.pdf",%]

ListLogPlot[{Table[{datapoints02[[k,1]],datapoints02[[k,2]]/datapoints02[[1,2]]},{k,1,Length[datapoints02]}],Table[{datapoints06[[k,1]],datapoints06[[k,2]]/datapoints06[[1,2]]},{k,1,Length[datapoints06]}],Table[{datapoints10[[k,1]],datapoints10[[k,2]]/datapoints10[[1,2]]},{k,1,Length[datapoints10]}]},PlotStyle->{Thick,Thick,Thick},PlotTheme->"Monochrome"(*,PlotLegends->{0.2,0.6,1.0}*),AxesLabel->{Style["K",20],Style[ToExpression["\\kappa_1",TeXForm,HoldForm],20]},TicksStyle->Directive[FontSize->18]]
Export["Fig4a.pdf",%]