Skip to main content

Posts

Showing posts from 2016

replacement - How to transform abstract finite group to permutation group?

It seems that Mathematica only has group functionality for permutation groups? Then there is a step to transform the abstract finite one to a permutation one. As an example, consider the following problem Suppose $G= $, try to compute $ab^2ab^2a$. It can proceed by viewing the elements of group just as a string, and the rule of group can be view as a rule of string replace. Any way, I would like to use NonCommutativeMultiply instead: Firstly, let us define the rule, basically, we only need to define the rule for $ab$ and $ba$, since whenever we know how to commute these two generators, then we can use the fact $a^2=1=b^3$ to simplify our expression. It is not hard to compute $ab=b^2a$, $ba=ab^2$: rel = {c___ ** a ** a ** d___ :> c ** i ** d, c___ ** b ** b ** b ** d___ :> c ** i ** d, c___ ** a ** b ** d___ :> c ** b ** b ** a ** d, c___ ** b ** a ** d___ :> c ** a ** b ** b ** d , a___ ** i :> a, i ** a___ :> a} Now, let us define the set of group elements

list manipulation - Emulating R data frame getters with UpValues

What's the best way to emulate R's data frames functionality? This includes the ability to select rows and columns in a 2-dimensional table by the string identifiers positioned typically in the first row or column, see this related question , (and thanks to Leonid Shifrin for cluing me to the term 'data frame.') Wrapping a 2D array in the symbol DF enables leveraging upvalues to essentially overload Part for the getters, here defined separately for singleton and list string arguments; matching uses Position. To select rows at level 1 of the array: Part[DF[expr_], r_String, rest___] ^:= expr[[Position[expr, r][[1, 1]], rest]] Part[DF[expr_], r_List /; And @@ StringQ /@ r, rest___] ^:= expr[[Position[expr, #][[1, 1]] & /@ r, rest]] Similarly to get columns at level 2 of the array: Part[DF[expr_], All, c_String] ^:= Transpose[expr][[Position[Transpose@expr, c][[1, 1]]]] Part[DF[expr_], All, c_List /; And @@ StringQ /@ c] ^:= Transpose[expr][[Position[Transpose@ex

memory - How to delete InterpolationFunction created by NDSolve/ParametricNDSolve?

I have faced a uncommon problem. The situation is that I need use ParametricNDSolve/ParametricNDSolveValue to solve ODEs system in the non-linear model fitting process. So every iteration generating a new parameter set and create a new InterpolatingFunction object. I guess the older InterpolatingFunction objects are not abandoned and it may be the reason for memory increase rapidly during the process. Now the problem is whether I can clear the objects at every iteration and how to do it. For a similar problem, I copy a 4D PDE example from MathGroup Archive : i = 0; sol = NDSolve[{D[u[t, x, y, z], t] == D[u[t, x, y, z], x, x] + D[u[t, x, y, z], y, y] + D[u[t, x, y, z], z, z], u[0, x, y, z] == 0, u[t, 0, y, z] == Sin[t], u[t, 40, y, z] == 0, u[t, x, 0, z] == Sin[t], u[t, x, 40, z] == 0, u[t, x, y, 0] == Sin[t], u[t, x, y, 40] == 0}, {u}, {t, 0, 100}, {x, 0, 40}, {y, 0, 40}, {z, 0, 40}, MaxSteps -> Infinity, MaxStepSize -> 1, EvaluationMonitor :

sorting - Import multiple files from folder in order of date added

I have a folder that contains certain files I want to import. In my case those are .png files. As all the files I want to import contain the phrase "accum" in their filenames, I can use SetDirectory[directory]; (*directory is string of the directory that contains the files*) importString = "./*" <> "accum" <> "*"; files = Import[importString]; to obtain a list consisting of all the images I wanted to import. Of course this list is sorted by filename. What can I do to tell Import that the files should not be sorted by "name" but by other options like "date added" (which is what I want)? Answer Likely something like the following: Import[#, "PNG"] & /@ SortBy[FileNames[importString], FileDate[#,"Creation"]&] although I suspect there might be a more terse approach.

equation solving - How to substitute the following conditions into an expression?

I have an expression: $p=a\;b\; x + b^2\; y + a\;c\; z$. I want to substitute $a\;b=1$, $b^2 = 2$ and $a\;c = 4$ to obtain $p = x + 2y + 4z$. How can I tell Mathematica to do that? I dont know how to start. Answer A reliable approach would use the third argument of Reduce as variables to eliminate (see Behavior of Reduce with variables as domain ) Reduce[{p == a b x + b^2 y + a c z, a b == 1, b^2 == 2, a c == 4}, {p}, {a, b, c}] p == x + 2 y + 4 z In the former editions of Mathematica ( ver <= 4 ) Reduce used the third argument for eliminating another variables. It can be still used this way though it is not documented anymore, however its trace could be found in SystemOptions["ReduceOptions"] . If Reduce didn't work this way one would exploit Solve (it is still supposed to eliminate variables ), e.g.: Apply[ Equal, Solve[{p == a b x + b^2 y + a c z, a b == 1, b^2 == 2, a c == 4}, {p}, {a, b, c}], {2}][[1, 1]] or Eliminate : Eliminate[{p == a b x + b^2 y

list manipulation - How do I sum the products of the elements of a vector?

I want to evaluate a sum like $\sum\limits_{j > i}^n q_i^k\, q_j^k$, where $n$ is the length of the vector $q^k$. The vector $q^k$ is, for example, $[q_1, q_2, q_3, q_4]$ and $i$ and $j$ correspond to the index of the vector elements. This was my attempt; I got errors :/ Answer vec = Table[q[i], {i, 4}] (* {q[1], q[2], q[3], q[4]} *) With[{l = Length[vec]}, Sum[vec[[i]] vec[[j]], {i, 1, l}, {j, i + 1, l}] ] (* q[1] q[2] + q[1] q[3] + q[2] q[3] + q[1] q[4] + q[2] q[4] + q[3] q[4] *) Plus @@ Times @@@ Subsets[vec, {2}] (* q[1] q[2] + q[1] q[3] + q[2] q[3] + q[1] q[4] + q[2] q[4] + q[3] q[4] *) Update For binary vectors (see OP's comment ), a fast approach would be: vec = RandomInteger[{0, 1}, 10^6]; Binomial[Total[vec], 2] // RepeatedTiming (* {0.0010, 124964252556} *) (* comparison of the result with Jim Baldwin's answer *) %[[2]] === (Total[vec^k]^2 - Total[vec^(2 k)])/2 (* True *)

Dynamic Delimiter in Manipulate

Consider the following Manipulate[a, Dynamic@Grid[{ {"Control 1", Control[{{a, 0, ""}, {1, 0}}]}, If[a == 1, {"Subcontrol 1", Control[{{aa, 0, ""}, {1, 0}}]}, Unevaluated[Sequence[]]], {"Control 2", Control[{{b, 0, ""}, {1, 0}}]} }, Spacings -> {Automatic, {2 -> 1}}, Dividers -> {False, {2 -> Manipulate`Dump`$delimitercolor}}, Alignment -> {{Right, Left}, Automatic}], ControlPlacement -> Left] My goal is to include a dynamic delimiter between the two controls. Since Delimiter doesn't seem to work within a Grid , an alternative is to use Spacings and Dividers with the right colouring, as previously discussed here . This solution, however, doesn't account for possible dynamic changes in the Control section. For example, if Control 1 is checked, we get As one can see, the divider doesn't change position. A solution would be to make the slight change Spacings ->

Hypergeometric function with a matrix argument

I am looking for the evaluation of a Hypergeometric function with a matrix argument as for example in Koev and Edelman or as showcased in this Wikipedia article . From what I understand from Mathematica 's documentation, it only accepts a scalar as the last argument. Answer Matrix functions in MMA First of all MMA does in general not support matrix arguments in its standard functions. Therefore there are special functions MatrixExp and MatrixPower available. But, as will be shown in this answer, it is possible to create user defined functions via infinite series, and it turns out that MMA is surprisingly well able in dealing with these matrix functions. The idea is simple: a function with a known series expansion can be transformed into a Matrix function letting z^k -> MatrixPower[z,k] . Hypergeometric matrix function in one variable Let us start with the well known function $2F1(a,b,c;z)$ which in MMA is Hypergeometric2F1[a,b,c,z] , and define matrix2F1[a_, b_, c_, z_] := S

Possible bug involving Dataset/Query and RightComposition

Bug introduced in 10.0.0 and fixed in 10.0.1 Short Version Why are the results of the following expressions different? Dataset[{1, 2}][f /* g, z] Dataset[<|x->1, y->2|>][Values /* f /* g, z] (* g[f[{z[1],z[2]}]] g[f[z[{1,2}]]] *) Long Version Consider the following dataset queries: Dataset[{1, 2}][f , z] Dataset[{1, 2}][f /* g, z] Dataset[{1, 2}][f /* g /* h, z] (* f[{z[1],z[2]}] g[f[{z[1],z[2]}]] h[g[f[{z[1],z[2]}]]] *) The results are as expected, with the ascending operators f , z , etc. being applied from the lowest level outward. But now consider these queries: Dataset[<|x->1, y->2|>][Values /* f, z] Dataset[<|x->1, y->2|>][Values /* f /* g, z] Dataset[<|x->1, y->2|>][Values /* f /* g /*h, z] (* f[{z[1],z[2]}] g[f[z[{1,2}]]] h[g[z[f[{1,2}]]]] *) This time, Values is used to ignore the keys in the assocation. The level 1 query operator is a composition of a descending operator with one or more ascending operators. The documentation

front end - Programming scripts to create and modify stylesheets: problems with contexts

I am trying to emulate the shift and shift + tab feature in the Outline.nb stylesheet that comes with Mathematica . I added the following code to a private stylesheet at the notebook level. Cell[StyleData["NUM"], CellDingbat->Cell[ TextData[{ CounterBox["NUM"], "."}]], CellMargins->{{80, 10}, {7, 7}}, ReturnCreatesNewCell->True, StyleKeyMapping->{"Tab" -> "SubNUM"}, CellGroupingRules->{"SectionGrouping", 50}, DefaultNewCellStyle->"SubNUM", DefaultReturnCreatedCellStyle->"NUM", ParagraphIndent->0, CounterIncrements->"NUM", CounterAssignments->{{"SubNUM", 0}}, MenuSortingValue->1200, MenuCommandKey->"1", FontFamily->"Times", FontColor->GrayLevel[0]] Cell[StyleData["SubNUM"], CellDingbat->Cell[ TextData[{ CounterBox["SubNUM"], "."}]], CellMargins->{{120, 10

graphics - Visualizing .xyz file in Mathematica

If i import a .xyz file in Mathematica that will show me the moelcule in the default way. But i want to visualise the .xyz file in the following manner. Answer Sticks and balls file = "ben.xyz"; atm = Import[file, "VertexTypes"]; xyz = Import[file]; coord = xyz[[1, 4, 1]]; balls = Cases[xyz[[1, 4, 2]], Sphere[x_, y_] -> {x, y}]; sticks= Cases[xyz[[1, 4, 2]], Cylinder[x_, y_] -> {x, y}]; Graphics[{EdgeForm[Black], FaceForm[White], {Thickness[#[[2]]/1000], Line[coord[[#[[1]]]][[All, 1 ;; 2]]]} & /@ sticks, {Disk[coord[[#[[1]]]][[1 ;; 2]], #[[2]]]} & /@ balls, Text[Style[atm[[#]], 20, Bold], coord[[#, 1 ;; 2]]] & /@ Range[Length[atm]]}]

Filtering of a Dataset based on a column name with underscore and a nested string pattern

I have the following dataset : ds=Dataset@{<|"rid" -> "#30:0", "objID" -> 1, "out_LTSnamAbbr" -> {"TMAP-TIR"}, "out_LTSnamPref" -> {"Topic Map - Term Information Resource"}, "out_LTSnamAlt" -> {"Index_term"}|>, <|"rid" -> "#30:1", "objID" -> 2, "out_LTSnamAbbr" -> {"TMAP-BIR"}, "out_LTSnamPref" -> {"Topic Map - Binary Information Resource"}, "out_LTSnamAlt" -> Missing["KeyAbsent", "out_LTSnamAlt"]|>, <|"rid" -> "#30:2", "objID" -> 3, "out_LTSnamAbbr" -> {"NAM-prefer"}, "out_LTSnamPref" -> {"Topic Map - Base name"}, "out_LTSnamAlt" -> Missing["KeyAbsent", "out_LTSnamAlt"]|>, <|"rid" -> "#30:3", "objID&qu

latex - How do I "start a reference system"?

When I attempt to insert a citation or specify a database to use as a bibliography, I get a message that I need to "start" the "reference system" "manually": I have no idea what a references system is, let alone how to start one ("manually" or otherwise). I've followed all of the steps in the documentation, to no avail. If it matters, I'm running OS 10.8.2, using a BibTeX file (which I've used without issues in other applications, e.g. LyX ) for my references, and I have a fill MacTeX installation in the standard locations on my machine. UPDATE: It matters . Answer At last (after 2.5 months), an answer from Wolfram "support" regarding my problems with this feature, documentation (and repeated insistence from Wolfram notwithstanding): It's only presently supported under Windows.

differential equations - Convert DE to first order?

I have the following set of second-order differential equations: eqs = {Derivative[2][x][t] + l0 Sin[ϕ[ t]] (Sin[θ[t]] (Derivative[1][θ][t]^2 + Derivative[1][ϕ][t]^2) - Cos[θ[t]] Derivative[2][θ][t]) - l0 Cos[ϕ[ t]] (2 Cos[θ[t]] Derivative[1][θ][ t] Derivative[1][ϕ][t] + Sin[θ[t]] Derivative[2][ϕ][t]) == 0, Derivative[2][y][t] + l0 Cos[θ[ t]] (-2 Sin[ϕ[t]] Derivative[1][θ][ t] Derivative[1][ϕ][t] + Cos[ϕ[t]] Derivative[2][θ][t]) - l0 Sin[θ[ t]] (Cos[ϕ[t]] (Derivative[1][θ][t]^2 + Derivative[1][ϕ][t]^2) + Sin[ϕ[t]] Derivative[2][ϕ][t]) == 0, l0 Cos[θ[t]] Derivative[1][θ][t]^2 + Derivative[2][z][t] + l0 Sin[θ[t]] Derivative[2][θ][t] == 0, l0 (-g m Cos[θ[t]] - l0 Cos[θ[t]] Sin[θ[t]] Derivative[1][ϕ][ t]^2 + Cos[θ[ t]] (-Sin[ϕ[t]] Derivative[2][x][t] + Cos[ϕ[t]] Derivative[2][y][t]) + Sin[θ[t]] Derivativ

polygons - Point lattice leading to triangle lattice

My main purpose is to eventually generate a triangle lattice from one original point, let's say the origin. So I want to start with the origin and generate 6 points around it, which are the vertices of an hexagon of side with length 1. Then, I would like to have those new six points to do the same and generate six other points (I know there will be repetitions). Then, I want those points to from a triangle lattice, but I am new to mathematica. So far I constructed a function that generates the vertices of an hexagon: h[x_, y_] := Point[Table[{Cos[2 Pi k/6] + x, Sin[2 Pi k/6] + y}, {k, 6}]] Now i am looking to generate the other poits from my new ones.

graphics3d - What are the possible ways of visualizing a 4D function in Mathematica?

I have a function $F$ that maps the xyz space to a set of reals, more clearly: $c = F[x,y,z]$ Where $c$,$x$,$y$ and $z$ are reals. What are the possible ways of visualizing this 3d function in Mathematica? (if possible, please post a how-to-do-it) Answer One possible way is to use Graphics3D with Point and color points by function value so it's like density plot 3d. For example, xyz = Flatten[ Table[{i, j, k}, {i, 1, 10, .35}, {j, 1, 10, .35}, {k, 1, 10, .35}], 2]; f[x_, y_, z_] := x^2 y Cos[z] Graphics3D[ Point[xyz, VertexColors -> (Hue /@ Rescale[f[##] & @@@ xyz])], Axes -> True, AxesLabel -> {x, y, z}] Another possible choice is just thinking one parameter as time variable and use Manipulate: Manipulate[Plot3D[f[x, y, z], {x, 1, 10}, {y, 1, 10}], {z, 1, 10}] There should be many other way to visualize 4d data, but it's really depending on what you want to see and how you want to visualize. Like amr suggested, you can also use Image3D

numerical value - Could I define 0^0 to be 1?

Many occasions, where I need to work numerically with functions. For a variable strictly between 0 and 1, during the optimization, it could become 0.^0 or 0^0, which then become indeterminate. Is there a way to define this 0^0=1? What are the possible down side of defining such relationship? Thanks! Answer With the help of @Michael E2 in my question Case $\frac{0}{0}$ In this case,you can define your function like this: func1[a_,b_]:=0 /;b==0 func1[a_,b_]:=a/b Test func1[0, 0] 1 Case$0.^0$ So you can use the /; to avoid $0.^0$ func2[x_,0]:=1/;x==0||x==0. func2[x_,y_:0]:=x^y Test func2[0, 0] 1 func2[0., 0] 1

calculus and analysis - Optimization over a bi-linear function

For the function $$ f(x_1,x_2;a_0,b_0)=\\\small\cases{\frac{1}{2}\left[x_1+x_2-x_1x_2+\left(\frac{-1+b_0(1-a_0)}{a_0}x_1+1\right)\left(\frac{-1+b_0(1-a_0)}{a_0}x_2+1\right)\right],\quad 0 \leq x_1\leq a_0,\, 0 \leq x_2\leq a_0\\ \frac{1}{2}\left[x_1+x_2-x_1x_2+\left(\frac{-1+b_0(1-a_0)}{a_0}x_1+1\right)b_0(-x_2+1)\right],\quad\quad\,\,\,\,\,\quad 0 \leq x_1\leq a_0,\, a_0 \leq x_2\leq 1\\ \frac{1}{2}\left[x_1+x_2-x_1x_2+b_0(-x_1+1)\left(\frac{-1+b_0(1-a_0)}{a_0}x_2+1\right)\right],\quad\quad\,\,\,\,\quad a_0 \leq x_1\leq 1,\, 0 \leq x_2\leq a_0\\ \frac{1}{2}\left[x_1+x_2-x_1x_2+b_0(-x_1+1)b_0(-x_2+1)\right],\quad\quad\quad\quad\quad\quad\quad a_0 \leq x_1\leq 1,\, a_0 \leq x_2\leq 1}$$ and $$f(x_1=x_2;a_0,b_0)=\cases{\frac{1}{2}\left[2x_1-{x_1}^2+\left(1+\frac{x_1(1+b_0(1-a_0))}{a_0}\right)^2\right],\quad 0 \leq x_1 \leq a_0\\\frac{1}{2}\left[2x_1+b_0^2(1-x_1)^2-x_1^2\right],\quad\quad\quad\quad\,\quad a_0 \leq x_1 \leq 1}$$ I would like to solve the following optimization problem: $$\

list manipulation - Problem with function inside brackets. Bug?

This code adds random digits to lists, and it works fine: a = {{1}, {2}, {3}}; Do[ j = RandomInteger[{1, Length[a]}]; AppendTo[a[[j]], RandomInteger[9]]; Print[a], {i, 5}]; (* {{1,7},{2},{3}} {{1,7},{2,2},{3}} {{1,7},{2,2},{3,9}} {{1,7},{2,2},{3,9,1}} {{1,7,2},{2,2},{3,9,1}} *) But if I replace the 'j' inside the [[]] with the definition j in the previous line, everything goes haywire: a = {{1}, {2}, {3}}; Do[ AppendTo[a[[RandomInteger[{1, Length[a]}]]], RandomInteger[9]]; Print[a], {i, 5}]; (* {{1},{1,7},{3}} {{1},{1,7},{1,4}} {{1,7},{1,7},{1,4}} {{1,7},{1,7},{1,7,9}} {{1,7,9},{1,7},{1,7,9}} *) Is this a bug or something I'm doing wrong? Answer Here's the issue. In the second (non-working) code, RandomInteger[{1, Length[a]}] is evaluated twice , as we can see by Trace ing the evaluation: SeedRandom[2] a = {{1}, {2}, {3}}; Trace[AppendTo[a[[RandomInteger[{1, Length[a]}]]], RandomInteger[9]], TraceInternal -> True] {RandomInteger[9

output formatting - Differently colored edges of same direction in a graph?

When I try to color edges in a graph that point to different vertices, it works properly: Graph[{Style[DirectedEdge[1, 2], Red], Style[DirectedEdge[2, 3], Blue]}] However, if I try to give different colors to two edges that involve the same vertices, it does not seem to work properly: Graph[{Style[DirectedEdge[1, 2], Red], Style[DirectedEdge[1, 2], Blue]}] We see that only one of the colors is registered. How should I apply the coloring properly in this case? Answer Modifying the answer in this Q/A : styles={Red, Blue, Green, Orange}; i=1; Graph[{a->b,a->b,a->b, a->b}, EdgeShapeFunction->({styles[[i++]],Arrow@#}&)] Update: dealing with more general cases Needs["GraphUtilities`"] styles1 = {Red, Blue, Green, Orange}; styles2 = { Purple, Cyan, Pink}; labels1 = {"A", "B", "C", "D"}; labels2 = { "E", "F", "G"}; Module[{i=1,j=1,i2,j2}, Graph[{a->b,a->b,a->b, a->b, a->c,c->

plotting - Using ListPointPlot3D to simulate 2D plots moving in time

I wrote a 1D solver for the heat equation $u_t=u_{xx}$, and I can animate the solution using normal ListPlot command, where the x-axis is the rod length, and the y-axis is $u(x,t)$. So for each time instance, I make a new 2D plot, etc... Another way to view this, is by 3D, where now each frame is a 2D plot at some time instance, and time travels down the page. This gives a nicer view of the diffusion of heat, I think. So, basically I have a number of 2D frames, and I want to display them in 3D, one frame after the other as time advances. Using MATLAB I use the mesh command with holdon , and here is a screen shot just to show what I mean from a school project I did: So time advances down the page as the simulation runs. I wrote the same thing using Mathematica for a demo, here is what I have to far, I got the general layout working using ListPointPlot3D , but I am not able to figure out how to tell Mathematica to make a curve on top of the points being displayed so I can get the same

function construction - Improve timings for calculating coeficients

Is there a quicker way of finding ther coefficients of something like $$(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}+x_{5}^{2}+x_{6}^{2}+x_{7}^{2}+x_{8}^{2}+x_{9}^{2}+x_{10}^{2})\times\\ (x_{1}+x_{2}+x_{3}+x_{4}+x_{5}+x_{6}+x_{7}+x_{8}+x_{9}+x_{10})^{8}$$ I have used s10 = Table[ToExpression["Subscript[x, " <> ToString[i] <> "]"], {i, 10}]; s10a = CoefficientRules@Expand[Plus @@ (s10^2)*(Plus @@ s10)^8]; Table[s10a[[n]], {n, Flatten[Position[s10a, #] & /@ (PadRight[#, 10, 0] & /@ IntegerPartitions@10), 1][[All, 1]]}] to find all of the coefficients of the distinct monomial terms using Mr.WIzard's suggestion below. Can I improve on this? Answer If interested in the evolution of this, read on, otherwise, skip to near the bottom for latest iteration... Give this a whirl, ran fine on a junky old netbook with 1 Gig that I keep at the local cigar shop: var = 10 pow = 12 result = With[{ar = Array[x, var]}, CoefficientArrays[Tr[ar^2]*Tr[ar]^po

expression manipulation - Count number of summands?

Consider fun = a+b+c+d/e If we want to get the number of summands, we can use Length[fun] , which properly gives 5 in this case. However, if fun contains only a single term fun = d/e Then applying Length[fun] gives 2 since now it actually counts the number of terms in the multiplication instead of summation. Therefore, Length is rather a hack than an actually reliable function to get the number of summands. Is there an efficient function that returns the number of summands reliably? Answer Best thing I can come up with is cntSummands[expr_]:=If[Head[expr]===Plus,Length[expr],If[expr === 0, 0, 1]] but this sounds like a terrible workaround. I am sure there are better ways?

list manipulation - Compile with array operations

When I use the function f with only the Module f[a1_, a2_, a3_]:= Module[.....] f[0,0,1200] works fine but gives errors when I use compile as follows. demand[n_, k_, qr_] := Min[k Vf, n capacity - qr]; supply[n_, k_, qr_] := Min[(n Kj - k) w, n capacity - qr]; Nsupply[n_, k_, qsum_, qrsum_, qr_] := Min[(n Kj - k) Vf - qsum - qrsum, n capacity - qr]; floExct[n_, ku_, k_, qsum_, qrsum_, qr_] := Min[demand[n, ku, qr], Nsupply[n, k, qsum, qrsum, qr]]; gamma[ku_, kd_] := Min[1, supply[L, kd, 0]/(demand[L, ku, 0] + 0.001)]; inflow[phi_, qin_] := (phi - \[Beta] qin) dx; density[qin_, qout_, qr_] := (qin - qout + qr)/Vf; Kj = 150.; w = 20.; Vf = 100.; capacity = 2500.; n = 30; m = 100; p = 27; RML = 18; \[Alpha] = 1200.; \[Beta] = 0.; L = 1.; delta = 1.; dt =4./3600.; dx = 1./9.; f= Compile[{{a1, _Integer}, {a2, _Integer}, {a3, _Integer}}, Module[{k0 = ConstantArray[0., n], Fk = Table[Table[0., {i1, n}], {j, 5}], kr = Table[Table[0., {i1, 1, p}], {i2, 1, n - 2}],