There are very few things I find unsatisfactory in L.C. Evans’ wonderful textbook on Partial Differential Equations; one of them is the illustration (on p.53 of the second edition) of the “heat ball”.

The heat ball is a region with respect to which an analogue of the mean value property of solutions to Laplace’s equation can be expressed, now for solutions of the heat equation. In the case of the Laplace’s equation, the regions are round balls. In the case of the heat equation, the regions are somewhat more complicated. They are defined by the expression

where is the fundamental solution of the heat equation

In the expressions above, the constant is the number of spatial dimensions; is the analogue of the radius of the ball, and in , the point is the center. Below is a better visualization of the heat balls: the curves shown are the boundaries in dimension , for radii between 0.75 and 4 in steps of 0.25 (in particular all the red curves have integer radii). In higher dimensions the shape is generally the same, though they appear more “squashed” in the direction.

1-dimensional heat balls centered at (0,5) for various radii. (Made using Desmos)

I’ve just done something that is, admittedly, rather silly. And I am somewhat surprised that it actually worked.

I managed to have TeXLive running on my new Samsung Galaxy Tab A, which runs Android.

It is painfully slow (speed is comparable to my old netbook from 2010). But in a pinch, I now know that it works. And that gives me some (minor) peace of mind.

Recently I started rethinking how I organize my incomplete and under development notes.

I know full well the inherent dangers of such an exercise in terms of my actual productivity. But now that I have completed my newest workflow I think I’ve finally found one that works well. (Fingers crossed.)

Before I describe what I do now, I’d like to document what I used to do, what I changed the last time, and why I am changing it again.

Joseph O’Rourke’s question at MathOverflow touched on an interesting characterization of geodesics in pseudo-Riemannian geometry, which was apparently originally due to Einstein, Infeld, and Hoffmann in their analysis of the geodesic hypothesis in general relativity. (One of my two undergraduate junior theses is on this topic, but I certainly did not appreciate this result as much when I was younger.) Sternberg’s book has a very good presentation on the theorem, but I want to try to give a slightly different interpretation in this post.

Geodesics and variation
One of the classical formulation of the criterion for a curve to be geodesic is that it is a stationary point of the length functional. Let be a Riemannian manifold, and let $latex: \gamma:[0,1]\to M$ be a mapping. Define the length functional to be

A geodesic then is a curve that is a critical point of under perturbations that fix the endpoints $\latex \gamma(0)$ and .

One minor annoyance about the length functional is that it is invariant under reparametrization of , and so it does not admit unique solutions. One way to work around this is to instead consider the energy functional (which also has the advantage of also being easily generalizable to pseudo-Riemannian manifolds)

It turns out that critical points of the energy functional are always critical points of the length functional. Furthermore, the energy functional has some added convexity: a curve is a critical point of the energy functional if it is a geodesic and that it has constant speed (in the sense that is independent of the parameter ).

The standard way to analyze the variation of is by first fixing a coordinate system . Writing the infinitesimal perturbation as , we can compute the first variation of :

Integrating the second term by parts we recover the familiar geodesic equation in local coordinates.

There is a second way to analyze the variation. Using the diffeomorphism invariance, we can imagine instead of varying while fixing the manifold, we can deform the manifold while fixing the curve . From the point of view of the energy functional the two should be indistinguishable. Consider the variation , which can be regarded as a vector field along which vanishes at the two end points. Let be a vector field on that extends . Then the infinitesimal variation of moving the curve in the direction should be reproducible by flowing the manifold by and pulling back the metric. To be more precise, let be the one parameter family of diffeomorphisms generated by the vector field , the first variation can be analogously represented as

By the definition of the Lie derivative we get the following characterizing condition for a geodesic:

Theorem
A curve is an affinely parametrized geodesic if and only if for every vector field vanishing near and , the integral

Noticing that , where is the Levi-Civita connection, we have that the above integral condition is equivalent to requiring

Using the boundary conditions and integrating by parts we see this also gives us, without passing through the local coordinate formulation, the geodesic equation

The Einstein-Infeld-Hoffmann theorem
The EIH theorem reads:

Theorem (EIH)
A curve is geodesic if and only if there exists a non-vanishing contravariant symmetric two tensor along such that for every vector field vanishing near and , the integral

(where is the induced length measure on ).

The EIH theorem follows immediately from the discussion in the previous section and the following lemma.

Lemma
A contravariant symmetric two tensor that satisfies the assumptions in the previous theorem must be proportional to .

Proof: Choose an orthonormal frame along for such that is tangent to . Write . Suppose . Then there exists a vector field such that and the symmetric part of is equal to . (We can construct by choosing a local coordinate system in a tubular neighborhood of such that . Then can be prescribed by its first order Taylor expansion in the normal direction to .) Let be a non-negative cut-off function and setting we note that since vanishes along . Therefore we have that the desired integral condition cannot hold. q.e.d.

References:

Shlomo Sternberg, Curvature in Mathematics and Physics

Einstein, Infeld, Hoffmann, “Gravitational Equations and the Problem of Motion”

A prominent academic, who happens not to be a mathematician, visited my home institution recently and gave a public address about the role of the university in the modern world. Most of what he said concerning our teaching mission are the usual platitudes about not being stuck in the past and making sure that our curricular content and learning objectives are aligned with what we would expect a 21st century college graduates to need.

It however bugged me to no end that the recurring example this particular individual returns to for something old-fashioned and “ought not be taught” is integration by parts; and he justifies this by mentioning that computer algebra systems (or even just google) can do the integrals faster and better than we humans can.

I don’t generally mind others cracking jokes at mathematicians’ expense. But this particular self-serving strawman uttered by so well-regarded an individual is, to those of us actually in the field teaching calculus to freshmen and sophomores, very damaging and disingenuous.

I happened to have just spent the entirety of last year rethinking how we can best teach calculus to the modern engineering majors. Believe me, students nowadays know perfectly well when we are just asking them to do busywork; they also know perfectly well that computer algebra systems are generally better at finding closed-form integral expressions than we can. Part of the challenge of the redesign that I am involved in is precisely to convince the students that calculus is worth learning in spite of computers. The difficulty is not in dearth of reason; on the contrary, there are many good reasons why a solid grounding of calculus is important to a modern engineering students. To give a few examples:

Taylor series are in fact important because of computers, since they provide a method of compactly encoding an entire function.

Newton’s method for root finding (and its application to, say, numerical optimization) is build on a solid understanding of differential calculus.

The entirety of the finite element method of numerical simulation, which underlies a lot of civil and mechanical engineering applications, are based on a variational formulation of differential equations that, guess what, only make sense when one understand integration by parts.

The notion of Fourier transform which is behind a lot of signal/image processing requires understanding how trigonometric functions behave under integration.

No, the difficulty for me and my collaborator is narrowing down a list of examples that we can not only reasonably explain to undergraduate students, but also have them have some hands-on experience working with.

When my collaborator and I were first plunged into this adventure of designing engineering-specific calculus material, one of the very first things that we did was to seek out inputs from our engineering colleagues. My original impulse was to cut some curricular content in order to give the students a chance to develop deeper understanding of fewer topics. To that end I selected some number of topics which I thought are old-fashioned, out-dated, and no longer used in this day and age. How wrong I was! Even something like “integration by partial fractions” which most practicing mathematicians will defer to a computer to do has its advocates (those who have to teach control theory insists that a lot of fundamental examples in their field can be reduced to evaluating integrals of rational functions, and a good grasp of how such integrals behave is key to developing a general sense of how control theory works).

In short, unlike some individuals will have you believe, math education is not obsolete because we all have calculators. In fact, I would argue the opposite: math education is especially pertinent now that we all have calculators. Long gone was the age where a superficial understanding of mathematics in terms of its rote computations is a valuable skill. A successful scientist or engineer needs to be able to effectively leverage the large toolbox that is available to her, and this requires a much deeper understanding of mathematics, one that goes beyond just the how but also the what and the why.

There are indeed much that can be done to better math education for the modern student. But one thing that shouldn’t be done is getting rid of integration by parts.

but with all the terms, whose decimal expansion includes the digit ‘9’, removed, in fact converges to some number below 80. The original proof is given in the Wikipedia article linked above, so I will not repeat it. But to make it easier to see the idea: let us first think about the case where the number is expressed in base 2. In base 2, all the positive integers has the leading binary bit being 1 (since it cannot be zero). Therefore there are no binary positive numbers without the bit ‘1’ in its expansion. So the corresponding series converges trivially to zero. How about the case of the bit ‘0’? The only binary numbers without any ‘0’ bits are

.

So the corresponding series actually becomes

So somewhere from the heavily divergent harmonic series, we pick up a rapidly converging geometric series. So what’s at work here? Among all the n-bit binary numbers, exactly 1 has all bits not being 0. So the density of these kinds of numbers decays rather quickly: in base 2, there are numbers that are exactly n-bit long. So if a number has a binary representation that is exactly n bits long (which means that ), the chances that it is one of the special type of numbers is . This probability we can treat then as a density: replacing the discrete sum by the integral (calculus students may recognize this as the germ of the “integral test”) and replacing the by the density , we get the estimate

.

Doing the same thing with the original Kempner series gives that the chances a n-digit number does not contain the digit nine to be

The length of the decimal expansion of a natural number is basically . So the density we are interested in becomes

From this we can do an integral estimate

The integral can be computed using that

to get

Notice that this estimate is much closer to the currently known value of roughly 22.92 than to the original upper bound of 80 computed by Kempner.

Kempner’s estimate is a heavy overestimate because he performed a summation replacing every n-digit long number that does not contain the digit 9 by ; this number can be many times (up to 9) times smaller than the original number. Our estimate is low because among the n-digit long numbers, the numbers that do not contain the digit 9 are not evenly distributed: they tend to crowd in the front rather than in the back (in fact, we do not allow them to crowd in the back because none of the numbers that start with the digit 9 is admissible). So if in the original question we had asked for numbers that do not contain the digit 1, then our computation will give an overestimate instead since these numbers tend to crowd to the back.

I use JabRef as my reference manager. In this post, however, I will discuss how we can abuse it to do some other things.

The problem

Let’s start with a concrete example: I keep a “lab notebook”. It is where I document all my miscellaneous thoughts and computations that come up during my research. Some of those are immediately useful and are collected into papers for publication. Some of those are not, and I prefer to keep them for future reference. These computations range over many different subjects. Now and then, I want to share with a collaborator or a student some subset of these notes. So I want a way to quickly search (by keywords/abstract) for relevant notes, and that compile them into one large LaTeX document.

Another concrete example: I am starting to collect a bunch of examples and exercises in analysis for use in my various classes. Again, I want to have them organized for easy search and retrieval, especially to make into exercise sheets.

The JabRef solution

The “correct” way to do this is probably with a database (or a document store), with each document tagged with a list of keywords. But that requires a bit more programming than I want to worry about at the moment.

JabRef, as it turns out, is sort of a metadata database: by defining a customized entry type you can use the BibTeX syntax as a proxy for JSON-style data. So for my lab notebook example, I define a custom type lnbentry in JabRef with

Required fields: year, month, day, title, file

Optional fields: keywords, abstract

I store each lab notebook entry as an individual TeX file, whose file system address is stored in the file field. The remaining metadata fields’ contents are self-evident.

(Technical note: in my case I actually store the metadata in the TeX file and have a script to parse the TeX files and update the bib database accordingly.)

For generating output, we can use JabRef’s convenient export filter support. In the simplest case we can create a custom export layout with the main layout file containing the single line

\\input{\file}

with appropriate begin and end incantations to make the output a fully-formed TeX file. Then one can simply select the entries to be exported, click on “Export”, and generate the appropriate TeX file on the fly.

(Technical note: JabRef can also be run without a GUI. So one can use this to perform searches through the database on the command line.)

One of my current research interests is in the geometry of constant-mean-curvature time-like submanifolds of Minkowski space. A special case of this is the study of cosmic strings under the Nambu-Goto action. Mathematically the classical (as opposed to quantum) behavior of such strings are quite well understood, by combining the works of theoretical physicists since the 80s (especially those of Jens Hoppe and collaborators) together with recent mathematical developments (such as the work of Nguyen and Tian and that of Jerrard, Novaga, and Orlandi). To get a better handle on what is going on with these results, and especially to provide some visual aids, I’ve tried to produce some simulations of the dynamics of closed cosmic strings. (This was also an opportunity for me to practice code-writing in Julia and learn a bit about the best practices for performance and optimization in that language.)

The code

After some false starts, here are some reasonably stable code.

function MC_corr3!(position::Array{Float64,2}, prev_vel::Array{Float64,2}, next_vel::Array{Float64,2}, result::Array{Float64,2}, dt::Float64)
# We will assume that data is stored in the format point(coord,number), so as a 3x1000 array or something.
num_points = size(position,2)
num_dims = size(position,1)
curr_vel = zeros(num_dims)
curr_vs = zeros(num_dims)
curr_ps = zeros(num_dims)
curr_pss = zeros(num_dims)
pred_vel = zeros(num_dims)
agreement = true
for col = 1:num_points #Outer loop is column
if col == 1
prev_col = num_points
next_col = 2
elseif col == num_points
prev_col = num_points - 1
next_col = 1
else
prev_col = col -1
next_col = col + 1
end
for row = 1:num_dims
curr_vel[row] = (next_vel[row,col] + prev_vel[row,col])/2
curr_vs[row] = (next_vel[row,next_col] + prev_vel[row,next_col] - next_vel[row,prev_col] - prev_vel[row,prev_col])/4
curr_ps[row] = (position[row,next_col] - position[row,prev_col])/2
curr_pss[row] = position[row,next_col] + position[row,prev_col] - 2*position[row,col]
end
beta = (1 + dot(curr_vel,curr_vel))^(1/2)
sigma = dot(curr_ps,curr_ps)
psvs = dot(curr_ps,curr_vs)
bvvs = dot(curr_vs,curr_vel) / (beta^2)
pssps = dot(curr_pss,curr_ps)
for row in 1:num_dims
result[row,col] = curr_pss[row] / (sigma * beta) - curr_ps[row] * pssps / (sigma^2 * beta) - curr_vel[row] * psvs / (sigma * beta) - curr_ps[row] * bvvs / (sigma * beta)
pred_vel[row] = prev_vel[row,col] + dt * result[row,col]
end
agreement = agreement && isapprox(next_vel[:,col], pred_vel, rtol=sqrt(eps(Float64)))
end
return agreement
end
function find_next_vel!(position::Array{Float64,2}, prev_vel::Array{Float64,2}, next_vel::Array{Float64,2}, dt::Float64; max_tries::Int64=50)
tries = 1
result = zeros(next_vel)
agreement = MC_corr3!(position,prev_vel,next_vel,result,dt)
for j in 1:size(next_vel,2), i in 1:size(next_vel,1)
next_vel[i,j] = prev_vel[i,j] + result[i,j]*dt
end
while !agreement && tries < max_tries
agreement = MC_corr3!(position,prev_vel,next_vel,result,dt)
for j in 1:size(next_vel,2), i in 1:size(next_vel,1)
next_vel[i,j] = prev_vel[i,j] + result[i,j]*dt
end
tries +=1
end
return tries, agreement
end

This first file does the heavy lifting of solving the evolution equation. The scheme is a semi-implicit finite difference scheme. The function MC_Corr3 takes as input the current position, the previous velocity, and the next velocity, and computes the correct current acceleration. The function find_next_vel iterates MC_Corr3 until the computed acceleration agrees (up to numerical errors) with the input previous and next velocities.

Or, in notations: MC_Corr3: ( x[t], v[t-1], v[t+1] ) --> Delta-v[t]
and find_next_vel iterates MC_Corr3 until Delta-v[t] == (v[t+1] - v[t-1]) / 2

The code in this file is also where the performance matters the most, and I spent quite some time experimenting with different algorithms to find one with most reasonable speed.

function make_ellipse(a::Float64,b::Float64, n::Int64, extra_dims::Int64=1) # a,b are relative lengths of x and y axes
s = linspace(0,2π * (n-1)/n, n)
if extra_dims == 0
return vcat(transpose(a*cos(s)), transpose(b*sin(s)))
elseif extra_dims > 0
return vcat(transpose(a*cos(s)), transpose(b*sin(s)), zeros(extra_dims,n))
else
error("extra_dims must be non-negative")
end
end
function perturb_data!(data::Array{Float64,2}, coeff::Vector{Float64}, num_modes::Int64)
# num_modes is the number of modes
# coeff are the relative sizes of the perturbations
numpts = size(data,2)
for j in 2:num_modes
rcoeff = rand(length(coeff),2)
for pt in 1:numpts
theta = 2j * π * pt / numpts
for d in 1:length(coeff)
data[d,pt] += ( (rcoeff[d,1] - 0.5) * cos(theta) + (rcoeff[d,2] - 0.5) * sin(theta)) * coeff[d] / j^2
end
end
end
nothing
end

This file just sets up the initial data. Note that in principle the number of ambient spatial dimensions is arbitrary.

using Plots
pyplot(size=(1920,1080), reuse=true)
function plot_data2D(filename_prefix::ASCIIString, filename_offset::Int64, titlestring::ASCIIString, data::Array{Float64,2}, additional_data...)
x_max = 1.5
y_max = 1.5
plot(transpose(data)[:,1], transpose(data)[:,2] , xlims=(-x_max,x_max), ylims=(-y_max,y_max), title=titlestring)
if length(additional_data) > 0
for i in 1:length(additional_data)
plot!(transpose(additional_data[i][1,:]), transpose(additional_data[i][2,:]))
end
end
png(filename_prefix*dec(filename_offset,5)*".png")
nothing
end
function plot_data3D(filename_prefix::ASCIIString, filename_offset::Int64, titlestring::ASCIIString, data::Array{Float64,2}, additional_data...)
x_max = 1.5
y_max = 1.5
z_max = 0.9
tdata = transpose(data)
plot(tdata[:,1], tdata[:,2],tdata[:,3], xlims=(-x_max,x_max), ylims=(-y_max,y_max),zlims=(-z_max,z_max), title=titlestring)
if length(additional_data) > 0
for i in 1:length(additional_data)
tdata = transpose(additional_data[i])
plot!(tdata[:,1], tdata[:,2], tdata[:,3])
end
end
png(filename_prefix*dec(filename_offset,5)*".png")
nothing
end

This file provides some wrapper commands for generating the plots.

include("InitialData3.jl")
include("MeanCurvature3.jl")
include("GraphCode3.jl")
num_pts = 3000
default_timestep = 0.01 / num_pts
max_time = 3
plot_every_ts = 1500
my_data = make_ellipse(1.0,1.0,num_pts,0)
perturb_data!(my_data, [1.0,1.0], 15)
this_vel = zeros(my_data)
next_vel = zeros(my_data)
for t = 0:floor(Int64,max_time / default_timestep)
num_tries, agreement = find_next_vel!(my_data, this_vel,next_vel,default_timestep)
if !agreement
warn("Time $(t*default_timestep): failed to converge when finding next_vel.")
warn("Dumping information:")
max_beta = 1.0
max_col = 1
for col in 1:size(my_data,2)
beta = (1 + dot(next_vel[:,col], next_vel[:,col]))^(1/2)
if beta > max_beta
max_beta = beta
max_col = col
end
end
warn(" Beta attains maximum at position $max_col")
warn(" Beta = $max_beta")
warn(" Position = ", my_data[:,max_col])
prevcol = max_col - 1
nextcol = max_col + 1
if max_col == 1
prevcol = size(my_data,2)
elseif max_col == size(my_data,2)
nextcol = 1
end
warn(" Deltas")
warn(" Left: ", my_data[:,max_col] - my_data[:,prevcol])
warn(" Right: ", my_data[:,nextcol] - my_data[:,max_col])
warn(" Previous velocity: ", this_vel[:,max_col])
warn(" Putative next velocity: ", next_vel[:,max_col])
warn("Quitting...")
break
end
for col in 1:size(my_data,2)
beta = (1 + dot(next_vel[:,col], next_vel[:,col]))^(1/2)
for row in 1:size(my_data,1)
my_data[row,col] += next_vel[row,col] * default_timestep / beta
this_vel[row,col] = next_vel[row,col]
end
if beta > 1e7
warn("time: ", t * default_timestep)
warn("Almost null... beta = ", beta)
warn("current position = ", my_data[:,col])
warn("current Deltas")
prevcol = col - 1
nextcol = col + 1
if col == 1
prevcol = size(my_data,2)
elseif col == size(my_data,2)
nextcol = 1
end
warn(" Left: ", my_data[:,col] - my_data[:,prevcol])
warn(" Right: ", my_data[:,nextcol] - my_data[:,col])
end
end
if t % plot_every_ts ==0
plot_data2D("3Dtest", div(t,plot_every_ts), @sprintf("elapsed: %0.4f",t*default_timestep), my_data, make_ellipse(cos(t*default_timestep), cos(t*default_timestep),100,0))
info("Frame $(t/plot_every_ts): used $num_tries tries.")
end
end

And finally the main file. Mostly it just ties the other files together to produce the plots using the simulation code; there are some diagnostics included for me to keep an eye on the output.

The results

First thing to do is to run a sanity check against explicit solutions. In rotational symmetry, the solution to the cosmic string equations can be found analytically. As you can see below the simulation closely replicates the explicit solution in this case.

The video ends when the simulation stopped. The simulation stopped because a singularity has formed; in this video the singularity can be seen as the collapse of the string to a single point.

Next we can play around with a more complicated initial configuration.

In this video the blue curve is the closed cosmic string, which starts out as a random perturbation of the circle with zero initial speed. The string contracts with acceleration determined by the Nambu-Goto action. The simulation ends when a singularity has formed. It is perhaps a bit hard to see directly where the singularity happened. The diagnostic messages, however, help in this regard. From it we know that the onset of singularity can be seen in the final frame:

The highlighted region is getting quite pointy. In fact, that is accompanied with the “corner” picking up infinite acceleration (in other words, experiencing an infinite force). The mathematical singularity corresponds to something unreasonable happening in the physics.

To make it easier to see the “speed” at which the curve is moving, the following videos show the string along with its “trail”. This first one again shows how a singularity can happen as the curve gets gradually more bent, eventually forming a corner.

This next one does a good job emphasizing the “wave” nature of the motion.

The closed cosmic strings behave like a elastic band. The string, overall, wants to contract to a point. Small undulations along the string however are propagated like traveling waves. Both of these tendencies can be seen quite clearly in the above video. That the numerical solver can solve “past” the singular point is a happy accident; while theoretically the solutions can in fact be analytically continued past the singular points, the renormalization process involved in this continuation is numerically unstable and we shouldn’t be able to see it on the computer most of the time.

The next video also emphasizes the wave nature of the motion. In addition to the traveling waves, pay attention to the bottom left of the video. Initially the string is almost straight there. This total lack of curvature is a stationary configuration for the string, and so initially there is absolutely no acceleration of that segment of the string. The curvature from the left and right of that segment slowly intrudes on the quiescent piece until the whole thing starts moving.

The last video for this post is a simulation when the ambient space is 3 dimensional. The motion of the string, as you can see, becomes somewhat more complicated. When the ambient space is 2 dimensional a point either accelerates or decelerates based on the local (signed) curvature of the string. But when the ambient space is 3 dimensional, the curvature is now a vector and this additional degree of freedom introduces complications into the behavior. For example, when the ambient space is 2 dimensional it is known that all closed cosmic strings become singular in finite time. But in 3 dimensions there are many closed cosmic strings that vibrate in place without every becoming singular. The video below is one that does however become singular. In addition to a fading trail to help visualize the speed of the curve, this plot also includes the shadows: projections of the curve onto the three coordinate planes.

Recently I have been playing around with the Julia programming language as a way to run some cheap simulations for some geometric analysis stuff that I am working on. So far the experience has been awesome.

A few random things …

Juno

Julia has a decent IDE in JunoLab, which is built on top of Atom. In terms of functionality it captures most of the sort of things I used to use with Spyder for python, so is very convenient.

Jupyter

Julia interfaces with Jupyter notebooks through the IJulia kernel. I am a fan of Jupyter (I will be using it with the MATLAB kernel for a class I am teaching this fall).

Plots.jl

For plotting, right now one of the most convenience ways is through Plots.jl, which is a plotting front-end the bridges between your code and various different backends that can be almost swapped in and out on the fly. The actual plotting is powered by things like matplotlib or plotlyJS, but for the most part you can ignore the backend. This drastically simplifies the production of visualizations. (At least compared to what I remembered for my previous simulations in python.)

Automatic Differentiation

I just learned very recently about automatic differentiation. At a cost in running time for my scripts, it can very much simplify the coding of the scripts. For example, we can have a black-box root finder using Newton iteration that does not require pre-computing the Jacobian by hand:

module NewtonIteration
using ForwardDiff
export RootFind
function RootFind(f, init_guess::Vector, accuracy::Float64, cache::ForwardDiffCache; method="Newton", max_iters=100, chunk_size=0)
### Takes input function f(x::Vector) → y::Vector of the same dimension and an initial guess init_guess. Apply Newton iteration to find solution of f(x) = 0. Stop when accuracy is better than prescribed, or when max_iters is reached, at which point a warning is raised.
### Setting chunk_size=0 deactivates chunking. But for large dimensional functions, chunk_size=5 or 10 improves performance drastically. Note that chunk_size must evenly divide the dimension of the input vector.
### Available methods are Newton or Chord
# First check if we are already within the accuracy bounds
error_term = f(init_guess)
if norm(error_term) < accuracy
info("Initial guess accurate.")
return init_guess
end
# Different solution methods
i = 1
current_guess = init_guess
if method=="Chord"
df = jacobian(f,current_guess,chunk_size=chunk_size)
while norm(error_term) >= accuracy && i <= max_iters
current_guess -= df \ error_term
error_term = f(current_guess)
i += 1
end
elseif method=="Newton"
jake = jacobian(f, ForwardDiff.AllResults, chunk_size=chunk_size, cache=cache)
df, lower_order = jake(init_guess)
while norm(value(lower_order)) >= accuracy && i <= max_iters
current_guess -= df \ value(lower_order)
df, lower_order = jake(current_guess)
i += 1
end
error_term = value(lower_order)
else
warn("Unknown method: ", method, ", returning initial guess.")
return init_guess
end
# Check if converged
if norm(error_term) >= accuracy
warn("Did not converge, check initial guess or try increasing max_iters (currently: ", max_iters, ").")
end
info("Used ", i, " iterations; remaining error=", norm(error_term))
return current_guess
end
end

This can then be wrapped in finite difference code for solving nonlinear PDEs!

I’ve just recently migrated to using NeoVim instead of traditional Vim. One of the nice features in NeoVim (or nvim) is that it now supports asynchronous job dispatch. This makes it a bit nicer to call external previewers for instance (otherwise the previewer may block the editing). So here are the latest LaTeX runtime code that I use, modified for NeoVim.

function Dvipreview()
let dviviewjob = jobstart(['xdvi', '-sourceposition', line(".")." ".expand("%"), expand("%:r") . ".dvi"])
endfunction
function PDFpreview()
let pdfviewjob = jobstart(['evince', expand("%:r") . ".pdf"])
endfunction
au BufRead *.tex call LaTeXStartup()
function LaTeXStartup()
set dictionary+=~/.config/nvim/custom/latextmp/labelsdictionary
set iskeyword=@,48-57,_,:
call SimpleTexFold()
set completefunc=CompleteBib
set completeopt=menuone,preview
runtime custom/latextmp/bibdictionary
call SetShortCuts()
endfunction
function SimpleTexFold()
exe "normal mz"
1
set foldmethod=manual
if search('\\begin{document}','nW')
1,/\\begin{document}/-1fold
if search('\\section','nW')
/\\section/1
endif
while search('\\section','nW')
.,/\\section/-1fold
/\\section/1
endwhile
.,$fold
endif
if search('\\begin{entry}','nW')
/\\begin{entry}/1
while search('\\begin{entry}','nW')
.,/\\begin{entry}/-1fold
/\\begin{entry}/1
endwhile
.,$fold
endif
exe "normal g`zzv"
endfunction
function SetShortCuts()
" Map <F2> to save and compile
imap <F2> ^[:w^M:! latex -src-specials % >/dev/null^M^Mi
" Map S-<F2> to save and compile as PDF
" apparently <S-F2> sends the same keycode as <F12>?
imap <F12> ^[:w^M:! pdflatex % >/dev/null^M^Mi
" Map <F3> to Dvipreview()
imap <F3> ^[:call Dvipreview()^M
" Map S-<F3> to PDFpreview()
" apparently <S-F3> = <F13>
imap <F13> ^[:call PDFpreview()^M
" Map <F4> to bibtex
imap <F4> ^[:! bibtex "%:r" >/dev/null^M^Mi
" Map <F5> to change the previous word into a latex \begin .. \end environment
imap <F5> ^[diwi\begin{^[pi<Right>}^M^M\end{^[pi<Right>}<Up>
" Map <F6> to 'escape the current \begin .. \end environment
imap <F6> ^[/\\end{.*}/e^Mi<Right>
" Map <F7> to search the labels dictionary for matching labels
imap <F7> ^[diwi\ref{^[pi<Right>^X^K
" Map <F8> to rebuild the labels dictionary
imap <F8> ^[:w^M:! ~/.config/nvim/custom/latexreadlabels.sh %^M^Mi
" Map <F9> to search using the bibs dictionary
imap <F9> ^[diwi\cite{^[pi<Right>^X^U
imap <S-Tab>C ^[diwi\mathcal{^[pi<Right>}
imap <S-Tab>B ^[diwi\mathbb{^[pi<Right>}
imap <S-Tab>F ^[diwi\mathfrak{^[pi<Right>}
imap <S-Tab>R ^[diwi\mathrm{^[pi<Right>}
imap <S-Tab>O ^[diwi\mathop{^[pi<Right>}
imap <S-Tab>= ^[diWi\bar{^[pi<Right>}
imap <S-Tab>. ^[diWi\dot{^[pi<Right>}
imap <S-Tab>" ^[diWi\ddot{^[pi<Right>}
imap <S-Tab>- ^[diWi\overline{^[pi<Right>}
imap <S-Tab>^ ^[diWi\widehat{^[pi<Right>}
imap <S-Tab>~ ^[diWi\widetilde{^[pi<Right>}
imap <S-Tab>_ ^[diWi\underline{^[pi<Right>}
endfunction

Pay attention that the control characters did not copy-paste entirely correctly in the SetShortCuts() routine. Those need to be replaced by the actual control-X sequences. The read labels shell script is simply

#!/bin/sh
grep '\label{' $1 | sed -r 's/.*\\label\{([^}]*)\}.*/\1/' > ~/.config/nvim/custom/latextmp/labelsdictionary

(I probably should observe the proper directory structure and dump the dictionary into ~/.local/share/ instead.)