text:
## NP-complete problems

In a previous post, I tried to explain the P vs NP problem. I gave a rough definition of the NP complexity class: NP is the class of the problems for which, for each “yes” instance, I can get convinced in polynomial time that the answer is yes by providing me with a proof that it is indeed the case. And, at the end of that post, I vaguely mentioned (without naming it) the notion of NP-complete problems.

The usual definition of a NP-complete problem is that it’s a problem which belongs to NP, and which is NP-hard. So now I need to explain what NP-hard means. Informally, it’s a problem that is at least as hard as all other problems of the NP class. Another informal way to say that is that it’s at least as hard as another NP-hard problem. The problem with that definition is the bootstrap problem: you’d need a first NP-hard problem to compare it to the others.

The archetypal NP-hard problem is SAT, which was the topic of my previous post. Since SAT is also in NP, SAT is also the archetypal NP-complete problem. It’s not a priori obvious that SAT is NP-hard. But one can prove that even 3-SAT is NP-hard, where 3-SAT represents the CNF instances of SAT where all clauses have exactly 3 literals. I’m not going to go through the proof, because that would imply that I’d explain a lot of things that I kind of don’t want to explain. For the people who are interested in the topic, it’s been proven in 1971, and it is known as the Cook-Levin theorem.

Now that we do have a NP-hard problem, we can find other ones by reduction. We reduce a NP-hard problem (for instance, 3-SAT), to another problem that we want to prove is NP-hard (call that problem A). To do that, we show that if we can solve A efficiently (that is to say, in polynomial time), then we can solve 3-SAT efficiently.

The idea of that kind of proof is to take an arbitrary instance of the NP-hard problem, and we transform it, in polynomial time, into an instance of A that allows to conclude on the original instance. Now suppose that we can solve problem A (that is to say, any instance of problem A) in polynomial time. Then we can take any instance of SAT, transform it into an instance of problem A, solve the instance of problem A, and get a conclusion. The transformation is done in polynomial time; the solving of A is in polynomial time; summing both yields again a polynomial, so we can solve any instance of SAT in polynomial time, so we can solve SAT in polynomial time. Easy.

NP-complete problems are, in a way, “complexity equivalent”: if we know how to solve one in polynomial time, then we can solve all of them in polynomial time as well, and we can solve all problems from NP in polynomial time (since NP-hard problems are at least as hard as any problem from NP). So, if we find a way to solve in polynomial time all the instances of any NP-complete problem… we proved that P = NP. And won a million dollars. And surprised a whole lot of people.

There is a large set of problems that are known to be NP-hard, or NP-complete if they are as well in NP. And there are people who look at… exotic problems, shall we say: let me give a few links for the people who are interested and may get a laugh out of it:

After this fun interlude, I’m going to give an example of reduction, so that you see that kind of proof works. I’m going to prove that the CLIQUE problem is NP-hard. It’s pretty much the basic example that everybody gives when they’re teaching that kind of things, but there must be a reason (I guess: it’s fairly easy to explain, and the reduction is also fairly easy compared to some others). What I’m explaining here is largely coming from the CLRS, that is to say the Cormen, Leiserson, Rivest and Stein, that is to say Introduction to Algorithms – that’s the reference book for A LOT of algorithms classes. (And I have a copy at home and a copy at work – you never know.)

The CLIQUE problem is described as follows: given a graph (a set of vertices and edges) $G$, does it contain a clique of size $k$? A clique of size $k$ is a set of $k$ vertices that are all connected to one another. For instance, this thing below is a clique of size 8. Also note that it does contain cliques of size 1 (single vertex), 2 (an edge), 3 (a triangle), 4, 5, 6, 7, since we can find sets of vertices of that size that are all connected to one another.

Let us reduce 3-SAT to that problem; since 3-SAT is NP-hard, CLIQUE will then be proven to be NP-hard as well. CLIQUE is in NP, because if I provide you with $k$ vertices, it’s possible to check that all these vertices are indeed connected to one another. So, if CLIQUE is NP-hard, CLIQUE is NP-complete.

To start with my reduction, I pick an arbitrary 3-SAT formula – more precisely, a 3-CNF-SAT formula (I explained the meaning of this in my post on SAT), which is a formula that has that look:

$(x_1 \vee x_2 \vee x_3) \wedge (\bar x_2 \vee x_4 \vee \bar x_5) \wedge (x_1 \vee \bar x_2 \vee \bar x_3) \wedge ...$

that is to say, a set of $m$ clauses connected by AND and composed of 3 literals connected by OR.

From there, we’re creating a graph that contains $3m$ vertices. Every vertex corresponds to a literal of the formula; vertices can be duplicated. For the above formula (truncated before the …), this yields the following vertices :

$x_1, x_2, x_3 ; \bar x_2, \bar x_4, \bar x_5 ; x_1, \bar x_2, \bar x_3$.

Then, we add edges almost everywhere. We don’t add edges in the “groups” that correspond to the clauses, and we also do not add edges between literals that are not compatible, that is to say inverse of each other. If I have two literals $x_1$ and $\bar x_1$, I’m not creating an edge between them. For the formula above, this is the graph in question:

And these are the edges that are NOT in the previous graph:

And now that we have that graph, we want to show that the SAT formula can be satisfied if and only if the graph (the first one) has a clique of size $m$, where $m$ is the number of clauses in the SAT formula. So I need to prove two things:

• if the formula can be satisfied, then the graph has a clique of size $m$,
• if the graph has a clique of size $m$, then the formula can be satisfied.

Let’s start with the first point. Suppose the formula can be satisfied. This means that we can assign a value to all the variables so that the formula is satisfied. This means that all clauses can be satisfied by this assignment. This also means that, for each clause, there’s a literal with value 1 (either a “positive” literal, for instance $x_1$ if the variable is assigned to 1, or a “negative” literal, for instance $\bar x_1$, if the variable is assigned to 0). Now remember that we created vertices grouped by “clause”, so for each clause, we can pick the vertex corresponding to that 1-valued literal (and if several literals have value 1 in a clause, we can pick an arbitrary one). Since we pick a vertex by clause, and we have $m$ clauses, we have $m$ vertices. We now want to prove that these $m$ vertices form a clique, which means that there is an edge between every pair of vertices of that set. We constructed the graph such that there is no edge between two vertices of a given clause, and we’re fine there, because we chose exactly one vertex per clause. Moreover, there is no edge between a literal and its negation – we’re also fine there, because we only picked literals that have value 1, and $x_i$ and $\bar x_i$ can’t both have value 1. These are the only conditions for which there is no edge between two vertices; which means that our $m$ vertices are all connected with each other, which yields a clique of size $m$, which is what we want.

Now, let’s look at the second point: if the graph has a clique of size $m$, then the formula can be satisfied. Suppose that the graph has a clique of size $m$. Since the vertices corresponding to the literals of a clause are not connected, this means that we have a vertex for each clause. We can give the value 1 to all the literals corresponding to these vertices. We cannot have a clique containing $x_i$ and $\bar x_i$, because there would be no edge between these two vertices, which goes against the definition of a clique, where all edges are present. So if a clique of size $m$ exists, this means that we found, for each clause, a literal whose value can be 1, and that all of these literals are compatible with each other. And so, we can satisfy the formula corresponding to the graph, which we also wanted to prove.

So, if we can solve CLIQUE in polynomial time, then we can solve 3-SAT in polynomial time; since 3-SAT is NP-hard, CLIQUE is NP-hard, so CLIQUE is NP-complete, which concludes the proof.

Showing that a problem is NP-complete is, in a way, an indicator that the problem in question is difficult, but this needs to be mitigated a lot. For one thing, it does not say anything about a specific instance of the problem. It does not say anything about a specific subset of instances either – let me explain that. If I say that CLIQUE is difficult, it doesn’t mean that, for example, deciding whether a triangle (a clique of size 3) is in a graph is difficult. I can take all sets of 3 vertices, look at whether they form a triangle, and conclude. There are approximately $n^3$ sets of 3 vertices in a graph with $n$ vertices (okay, there’s exactly $\displaystyle \frac{n(n-1)(n-2)}{6}$ – anyway, roughly $n^3$), so I can actually decide that in polynomial time (because I’m doing $n^3$ operations which are basically checking if three edges are in a graph). So I can decide 3-CLIQUE in polynomial time. Well, I’m not going to be a millionaire with that, because the CLIQUE problem is wider than just 3. I can also decide 1000-CLIQUE (clique of size 1000) in polynomial time with the same principle. Well, it’s a polynomial of degree 1000, but who cares 😛

But, in the general case, I cannot decide whether a graph over $n$ vertices contains a clique of $\displaystyle \frac{n}{2}$, or $\displaystyle \frac{n}{10000}$ vertices, or even $\log n$ vertices in polynomial time with this algorithm that looks at all groups of size $k$ (so, here, $k$ would be equal to $\displaystyle \frac{n}{2}$, $\displaystyle \frac{n}{10000}$, or $\log n$) and that looks at the edges of the group, because that would mean doing roughly $n^{n/2}$, $n^{n/10000}$ and $n^{\log n}$ operations, respectively, and that in any case it’s not polynomial. And, generally speaking, nobody knows if it’s possible to do that in polynomial time, since CLIQUE is NP-complete and we don’t know if P is equal to NP. Many people think that not, but nobody knows.

Even that does not say anything about the difficulty of a specific instance. If I’m given a graph that contains $\displaystyle \frac{n(n-1)}{2}$ edges for $n$ vertices, then I can decide pretty easily that it does contain a clique of size $\displaystyle \frac{n}{2}$, $\displaystyle \frac{n}{10000}$ and $\log n$. That’s because a graph that has $\displaystyle \frac{n(n-1)}{2}$ for $n$ vertices contains all the possible edges for the graph: it’s a clique on $n$ vertices… which means that every subgraph over $k$ vertices, for all $k \leq n$, is a clique of size $k$. Same thing if I’m given a graph over $n$ vertices with 0 edge, it’s not going to be hard to say that it doesn’t contain a clique of size $\displaystyle \frac n 2$.

SAT is also NP-hard, but can solve 2-CNF-SAT instances (where all clauses contain 2 literals) in polynomial, and even linear time. There are also SAT instances that are deemed “easy”, for instance the ones that have only few dependencies between the clauses: we know that we can to solve them, and we even know how (for the people who are interested in what I mean exactly by what I mean by “a few dependencies”, I’m thinking specifically about “the ones that are in the scope of the Lovász Local Lemma“, whose application to SAT is not trivial from the Wikipedia definition, but which might be clearer here).

Generally speaking, it may even be pretty difficult to create NP-hard instances of problems that we cannot solve in polynomial time. By the way, it’s not easy to show that we don’t know how to solve a given instance in polynomial time. We can find instances problem instances that are hard (i.e. require exponential time) for a given algorithm, and easy for another algorithm…

And that’s a bit of what happens in real life: we have instance of problems for which we know that the general case is NP-hard (or worse), and that we still solve without great difficulty. It can be that the problem is “small” (it’s not a huge problem to put an exponential algorithm on a problem of size 4), but it can also be that the problem belongs to an “easy” subset of instances. Considering that an instance of a problem is hard because the general problem is hard, and giving up on it (because we won’t know how to do it) may not be a good idea – it may be that we forgot to consider an easy restriction of the problem!

As usual, “it’s not that simple”. I’ve spent a fair amount of time with theory, and I tend to look at the problems themselves more than at given instances. But it’s sometimes nice to remember that stuff that’s “theoretically hard” can end up being “practically easy” – it’s a nice change from the things that are theoretically easy but that we don’t know how to do practically… 🙂


text:
# inter arrival time ( continuous exponential distribution)

The customers arrive at store with exponentially distributed interarrival times with a mean of 30 minutes. What is the probability that interarrival time between two successive customers is 2 hours or more?

My attempt: 2 customers arrive every hour, therefore λ=2.

P(T≥2)=∫ λ^(-λt) dt ( where $λ=2$, upperbound is infinity, and lowerbound is 2).This gives an answer of=$0.018$

however the answer in the answer key is $0.18$ can someone tell me what I did wrong. Thank you very much

The CDF of exponential distribution is $F_T(t)=1-e^{-\lambda t}$. If described in the unit of hours, we have $\lambda^{-1}=0.5$ and $\lambda=2$. The question asks for $\Pr(T>2)$ which is given by $1-F_T(2)$ or $e^{(-2)(2)}=0.018$, and your answer is correct.

format: plain


text:
# What is the total enthalpy for "C"(g) -> "C"^(4+)(g) + 4e^(-)?

Jun 6, 2017

$\text{148.0262 eV}$, or $\text{14280.78 kJ/mol}$.

Well, it's very unusual, but I guess we can do it... Break it into four steps.

The ionization energy is the energy put into the gaseous atom to remove an electron.

${\text{C"(g) -> "C}}^{+} \left(g\right) + {e}^{-}$, $\Delta {H}_{I {E}_{1}} = \text{11.260300 eV}$

${\text{C"^(+)(g) -> "C}}^{2 +} \left(g\right) + {e}^{-}$, $\Delta {H}_{I {E}_{2}} = \text{24.3846 eV}$

${\text{C"^(2+)(g) -> "C}}^{3 +} \left(g\right) + {e}^{-}$, $\Delta {H}_{I {E}_{3}} = \text{47.88778 eV}$

${\text{C"^(3+)(g) -> "C}}^{4 +} \left(g\right) + {e}^{-}$, $\Delta {H}_{I {E}_{4}} = \text{64.49351 eV}$

where the ionization energies were taken from NIST (most reliable for light elements).

(The electron volt, $\text{eV}$ is a unit of energy needed to accelerate the electron through a $\text{1 V}$ potential difference.)

The change in enthalpy, $\Delta H$, is additive, and a state function (a function with an initial and final state independent of the path).

So, we can build a series of steps that each have their own $\Delta H$, without really caring what particular path we take, as long as the initial and final states are what we want.

$\Delta {H}_{C \left(g\right) \to {C}^{4 +} \left(g\right)} = \Delta {H}_{I {E}_{1}} + \Delta {H}_{I {E}_{2}} + \Delta {H}_{I {E}_{3}} + \Delta {H}_{I {E}_{4}}$

$= \text{11.260300 eV" + "24.3846 eV" + "47.88778 eV" + "64.49351 eV}$

$=$ $\text{148.0262 eV}$

which is actually quite large, but that makes sense, right? You're ionizing carbon FOUR times, which is probably nearly impossible to accomplish in real life... (Even if it was, it would be short-lived.)

If you don't see how large this is, let's try it in $\text{kJ/mol}$.

color(blue)(DeltaH_(C(g) -> C^(4+)(g))) = 148.0262 cancel"eV" xx (1.602 xx 10^(-19) cancel"J")/cancel"eV" xx "1 kJ"/(1000 cancel"J") xx 6.0221413 xx 10^(23) ("mol C")^(-1)

$=$ $\textcolor{b l u e}{\text{14280.78 kJ/mol}}$

Compare that to the following ionization energies:

color(white)([(" ", color(black)("kJ/mol")),(color(black)("H"), color(black)(1312.0)),(color(black)("Na(2nd)"),color(black)(4562.4)),(color(black)("Mg(3rd)"),color(black)(7732.6)),(color(black)("Al(4th)"),color(black)(11577))])

It exceeds the fourth ionization energy of aluminum, which is the removal of its $2 s$ core electron, by $\text{2704 kJ/mol}$, which is also about $27$ times the enthalpy of the average chemical bond.


text:
# Question #9cdcc

Oct 15, 2017

$4 x$

#### Explanation:

well, the radical signs are unnecessary. $\sqrt{1} = 1$

...so you can rewrite as:

$1 - {x}^{2} - x - {x}^{2} + x$

...and I think this simplifies to:

$1 - 2 {x}^{2}$

And the derivative of this would b:

$- 4 x$

GOOD LUCK

Oct 15, 2017

Given: $\sqrt{1} - {x}^{2} - \frac{x}{\sqrt{1}} - {x}^{2} + x$

Use the substitution, $\sqrt{1} = 1$:

$1 - {x}^{2} - x - {x}^{2} + x$

Combine like terms:

$1 - 2 {x}^{2}$

Differentiate:

$\frac{d \left(1 - 2 {x}^{2}\right)}{\mathrm{dx}} = \frac{d \left(1\right)}{\mathrm{dx}} - \frac{d \left(2 {x}^{2}\right)}{\mathrm{dx}}$

The first term is 0 because the derivative of a constant is 0:

$\frac{d \left(1 - 2 {x}^{2}\right)}{\mathrm{dx}} = - \frac{d \left(2 {x}^{2}\right)}{\mathrm{dx}}$

Use the linear property of the derivative:

$\frac{d \left(1 - 2 {x}^{2}\right)}{\mathrm{dx}} = - 2 \frac{d \left({x}^{2}\right)}{\mathrm{dx}}$

Use the power rule, $\frac{d \left({x}^{n}\right)}{\mathrm{dx}} = n {x}^{n - 1}$:

$\frac{d \left(1 - 2 {x}^{2}\right)}{\mathrm{dx}} = - 4 x$


text:
up next prev ptail tail

## Article 9Variations of the Hyperbolic Functions.

Since the values of the hyperbolic functions depend only on the sectorial measure, it is convenient, in tracing their variations, to consider only sectors of one half of a rectangular hyperbola, whose conjugate radii are equal, and to take the principal axis $OA$ as the common initial line of all the sectors. The sectorial measure $u$ assumes every value from $-\infty$, through $0$, to $+\infty$, as the terminal point $P$ comes in from inﬁnity on the lower branch, and passes to inﬁnity on the upper branch; that is, as the terminal line $OP$ swings from the lower asymptotic position $y=-x$, to the upper one, $y=x$. It is here assumed, but is proved in Art. 17, that the sector $AOP$ becomes inﬁnite as $P$ passes to inﬁnity.

Since the functions $coshu,sinhu,tanhu$, for any position of $OP$, are equal to the ratios of $x,y,t$, to the principal radius $a$, it is evident from the ﬁgure that

 $cosh0=1,\phantom{\rule{1em}{0ex}}sinh0=0,\phantom{\rule{1em}{0ex}}tanh0=0,$ (15)

and that as $u$ increases towards positive inﬁnity, $coshu,sinhu$ are positive and become inﬁnite, but $tanhu$ approaches unity as a limit; thus

 $cosh\infty =\infty ,\phantom{\rule{1em}{0ex}}sinh\infty =\infty ,\phantom{\rule{1em}{0ex}}tanh\infty =1.$ (16)

Again, as $u$ changes from zero towards the negative side, $coshu$ is positive and increases from unity to inﬁnity, but $sinhu$ is negative and increases numerically from zero to a negative inﬁnite, and $tanhu$ is also negative and increases numerically from zero to negative unity; hence

 (17)

For intermediate values of $u$ the numerical values of these functions can be found from the formulas of Arts. 16, 17, and are tabulated at the end of this chapter. A general idea of their manner of variation can be obtained from the curves in Art. 25, in which the sectorial measure $u$ is represented by the abscissa, and the values of the functions $coshu$, $sinhu$, etc., are represented by the ordinate.

The relations between the functions of $-u$ and of $u$ are evident from the deﬁnitions, as indicated above, and in Art. 8. Thus

 (18)
Prob. 12.
Trace the changes in $sechu,cothu,cschu$, as $u$ passes from $-\infty$ to $+\infty$. Show that $sinhu,coshu$ are inﬁnites of the same order when $u$ is inﬁnite. (It will appear in Art. 17 that $sinhu,coshu$ are inﬁnites of an order inﬁnitely higher than the order of $u$.)
Prob. 13.
Applying eq. (12) to ﬁgure, page §, prove $tanh{u}_{1}=tanAOP$.

 up next prev ptail top


text:
1. ## Analyticity of z

I was reading through my text, and it stated that $\displaystyle z^n$ (for some integer n) was not analytic at n = 0. Why is this the case? Usually I'd substitute z = x + iy and then use the Cauchy-Riemann equations to check, but (x + iy)^n makes it difficult to separate the real and imaginary parts.

Thanks.

2. Originally Posted by Glitch
I was reading through my text, and it stated that $\displaystyle z^n$ (for some integer n) was not analytic at n = 0. Why is this the case? Usually I'd substitute z = x + iy and then use the Cauchy-Riemann equations to check, but (x + iy)^n makes it difficult to separate the real and imaginary parts.

Thanks.
Do you mean when $\displaystyle z=0$ not n?

3. Sorry, I should probably provide some context.

The question I got this from is the following:

Let C be the circle of radius R centre the origin, taken anticlockwise. Find $\displaystyle \int_C{z^n \dz}$ for integer values of n.

4. Originally Posted by Glitch
Sorry, I should probably provide some context.

The question I got this from is the following:

Let C be the circle of radius R centre the origin, taken anticlockwise. Find $\displaystyle \int_C{z^n \dz}$ for integer values of n.
Well you can parametrize the circle at the origin by

$\displaystyle z=Re^{i\theta},\theta \in [0,2\pi] \implies dz=Rie^{i\theta}d\theta$

When

$\displaystyle n=0 \implies z^n=1$ so you get the integral

$\displaystyle \oint dz=Ri\int_{0}^{2\pi}e^{i\theta}d\theta=0$

In fact you can use this for every integer n but be careful with the case n=-1

5. I just don't understand why it is not analytic at n = 0 (at least, according to this text).

6. Originally Posted by Glitch
I just don't understand why it is not analytic at n = 0 (at least, according to this text).
It should be when n=0 we get

$\displaystyle z^0=\begin{cases}1, z \ne 0 \\ \text{undefined},z=0 \end{cases}$

Technically $\displaystyle 0^0$ is undefined, but if we check the limit along any path of

$\displaystyle \lim_{z \to 0}z^0=1$ So by the Riemann removable singularity theorem the function can be analytically extended to all of the complex plane

Removable singularity - Wikipedia, the free encyclopedia

Note that holomorphic function is a function that is analytic on all of the complex plane.


text:
# An inequality on $C^1$ periodic functions

Suppose $f \in C^1(\mathbb{R})$ and $f(x + 1) = f(x) \ \forall x \in \mathbb{R}$. Show that $$||f||_{\infty} \leq \int_0^1|f| + \int_0^1|f'|.$$

I have tried using techniques in Fourier Analysis such as Parseval's Identity and expanding $f$ since $f$ is periodic and $C^1$; I have also tried Mean Value Theorem, but I have yet to produce anything useful. Both answers and hints are appreciated and thanks in advance!

• What does the mean value theorem tell you about the first integral? How does the second integral relate to the distance between $f(x)$ and $f(y)$ for $x,\,y \in [0,\,1]$? – Daniel Fischer Jul 8 '13 at 14:15

1. $\|f\|_{\infty}=\max_{0\le x\le 1}|f(x)|=|f(x_0)|$ for some $x_0\in[0,1]$
2. $|f(y_0)|=\min_{0\le y\le 1}|f(y)|$ for some $y_0\in[0,1]$(Note $|f(y_0)|\le \int_0^1|f(x)|\;\mathbb{d}x$)
3. $|f(x_0)|=|f(x_0)-f(y_0)+f(y_0)|=|\int_{y_0}^{x_0}f^{\prime}(x)\;\mathbb{d}x+f(y_0)|\le \int_0^1|f^{\prime}(x)|\;\mathbb{d}x+|f(y_0)| \\ \le \int_0^1|f^{\prime}(x)|\;\mathbb{d}x+\int_0^1|f(x)|\;\mathbb{d}x$
By the mean-value theorem there is a $0 \leq \xi \leq 1$ such that $$\int_{0}^{1} f(x) dx = f(\xi)$$ Therefore $$f(x) = f(x) - f(\xi) + \int_{0}^{1} f(y) dy = \int_{\xi}^{x} f'(y) dy + \int_{0}^{1} f(y) dy$$ Taking absolute values we get $$|f(x)| \leq \int_{0}^{1} |f'(y)| dy + \int_{0}^{1} |f(y)| dy$$ This inequality can be refined. I think you can put an $\tfrac 12$ in front of the derivative term.


text:
# DETERMINATION OF A TWO-DIMENSIONAL HEAT SOURCE : UNIQUENESS, REGULARIZATION AND ERROR ESTIMATE

Abstract : Let $Q$ be a heat conduction body and let $\varphi = \varphi(t)$ be given. We consider the problem of finding a two-dimensional heat source having the form $\varphi(t)f(x,y)$ in $Q$. The problem is ill-posed. Assuming $\partial Q$ is insulated and $\varphi \not\equiv 0$, we show that the heat source is defined uniquely by the temperature history on $\partial Q$ and the temperature distribution in $Q$ at the initial time $t = 0$ and at the final time $t = 1$. Using the method of truncated integration and the Fourier transform, we construct regularized solutions and derive explicitly error estimate.
Keywords :
Type de document :
Article dans une revue
Journal of Computational and Applied Mathematics, Elsevier, 2006, 191, pp.50-67

https://hal.archives-ouvertes.fr/hal-00009236
Contributeur : Alain Pham Ngoc Dinh <>
Soumis le : jeudi 29 septembre 2005 - 16:09:10
Dernière modification le : jeudi 3 mai 2018 - 15:32:06
Document(s) archivé(s) le : jeudi 1 avril 2010 - 22:35:02

### Identifiants

• HAL Id : hal-00009236, version 1

### Citation

Dang Duc Trong, Pham Hoang Quan, Alain Pham Ngoc Dinh. DETERMINATION OF A TWO-DIMENSIONAL HEAT SOURCE : UNIQUENESS, REGULARIZATION AND ERROR ESTIMATE. Journal of Computational and Applied Mathematics, Elsevier, 2006, 191, pp.50-67. 〈hal-00009236〉

### Métriques

Consultations de la notice

## 247

Téléchargements de fichiers


text:
summation related issues & queries in MathXchanger

Relabelling indices in multiple summations

Updated November 13, 2018 19:20 PM

Sum of $\lfloor k^{1/3} \rfloor$

Updated November 13, 2018 09:20 AM

Interesting proofs of this sum

Updated November 13, 2018 08:20 AM

Need to solve a sum

Updated November 10, 2018 16:20 PM

How to prove equality

Updated November 10, 2018 12:20 PM

Confusion with Maclaurin polynomials

Updated November 10, 2018 05:20 AM

Finding sum of the infinite series.

Updated November 09, 2018 17:20 PM

Campbell Theorem for n-th power

Updated November 08, 2018 21:20 PM

Stuck at summation of exponential.

Updated November 08, 2018 15:20 PM

nth power of Summation

Updated November 08, 2018 13:20 PM

Concrete Mathematics: Quicksort analysis

Updated November 08, 2018 01:20 AM

Prove a summation equals 1

Updated November 07, 2018 22:20 PM

Double sum identity

Updated November 07, 2018 16:20 PM

Optimizing a sum function for computing

Updated November 07, 2018 15:20 PM

Help to simplify a summation

Updated November 07, 2018 05:20 AM

Finding sum of a finite series.

Updated November 05, 2018 19:20 PM

sum $\vert i-j\vert^{\alpha}$

Updated November 05, 2018 14:20 PM

Showing Page 1 of 0


text:
# Approximating Geodesics in a half edge DS, how can I refine my mesh to get good approximations

I implemented Djikstra's shortest path algorithm to approximate Geodesics on arbitrary meshes. Djikstra's works, but I noticed a problem inherent to the discretization of my meshes.

Consider the following figure squence:

...

This is my current refinement algorithm which is the easiest/standard face subdivision. Now consider the approximation of a geodesic in 2 points:

The blue point is where I think the actual geodesic intersects that edge, which is quite far from where the approximated geodesic passes. However that path ISN'T wrong.

Consider a square grid. The distance between any 2 points in the grid is the manhattan distance |x| + |y|.

So as far as Djikstra's is concerned, a path that goes all the way down and then to the left has the same length as a path that goes diagonally in a staircase pattern. Refining the mesh won't change the distance either. In other words, the limit of the shortest path found by Djikstra in a regular square grid as the size of the squares goes to 0 is NOT the straigtht line connecting 2 points.

Now the actual question, does anyone know of a way to subdivide my surface that is fairly straightforward but will actually converge to the geodesic?

• Intriguing problem! If I am not mistaken, the situation is even worse than you expected in the sense that any path going through the blue point should be a little bit longer than the path you found with Dijkstra's algorithm. – user9485 Jul 21 at 14:15
• The analytic path should be shorter, but the approximation might be longer – Makogan Jul 21 at 15:45
• Yes, this is what I meant to say. I think you found the (unique) shortest path between the two red points, among all paths that move along the triangles' edges. – user9485 Jul 21 at 15:54
• In this specific case there's just one, but in general there might be multiple, including infinite. – Makogan Jul 21 at 17:02
• How can there be infinitely many shortest paths? – user9485 Jul 21 at 18:30

## 1 Answer

As you pointed out, the problem here is the discretization/subdivision of the mesh. If your mesh was made of quadrilaterals instead of triangles, the obvious subdivision strategy would be to split each quadliteral into four equally sized smaller quadliterals:

$$\hspace{2cm}$$ $$\hspace{2cm}$$

For any two points $$P_1$$ and $$P_2$$, Dijkstra's algorithm would yield sets of shortest paths between these points $$P_1$$ and $$P_2$$. The more you refine the discretization with this subdivision strategy, the more shortest paths you'll find between the two points. However, intuitively it is clear that for every subdivision level $$l\in\mathbb{N}$$ you could pick one of these shortest paths $$p_l$$ such that the sequence $$(p_l)_{l\in\mathbb{N}}$$ converges to the actual geodesic between $$P_1$$ and $$P_2$$ (with respect to some norm that has to be specified, for example the supremum norm).

Unfortunately, the same is not true for the standard subdivision strategy of splitting a triangle into four smaller triangles, which you have proven with your example. I believe that, at its core, the problem is that there is no way to reach the center of the triangle with a straight line from each of its edges. This can be achieved by splitting a triangle in each subdivision step into 6 smaller triangles like this:

$$\hspace{2cm}$$ $$\hspace{2cm}$$

I do not have a proof that this subdivision is more useful to compute geodesics with Dijkstra's algorithm, but it seems quite likely to me. I would be very interested to see what your results look like with this subdivision strategy! However, no matter what you do, in the end you may or may not end up with sets of shortest paths instead of a single shortest path. In this case you will need some kind of heuristic or additional algorithm to decide which path resembles the true geodesic most closely.

• I think you have it wrong, given the square subdivision you gave, a refinement of that will not converge towards the geodesic, because the distance between 2 points remains constant. I explained that on my original question : p – Makogan Jul 21 at 17:03
• Your triangle scheme has given me and Idea, so I will test that out and get back to you – Makogan Jul 21 at 17:04
• This is why I wrote sets of shortest paths. You always need to make sure you calculate all shortest paths and not just one random path. Then you indeed get convergence as I explained above. – user9485 Jul 21 at 18:40
• I am going to have to ask another question here before I get back to you – Makogan Jul 22 at 4:22
• This is not what I ended up doing but it;s probably better, my solution was just to start randomly splitting edges (I just did it this way because it;s easier in my current implementation) but this will probably work too. – Makogan Jul 22 at 4:46


text:
# Conditional PDF question -- I think anyway...

Tags:
1. Sep 27, 2015

### whitejac

1. The problem statement, all variables and given/known data
Suppose you take a pass-fail test repeatedly. Let Sk be the event that you are successful in your kth try, and Fk be the event that you fail the test in your kth try. On your first try, you have a 50% chance of passing the test.

P(S1)=1−P(F1)=1/2.
Assume that as you take the test more often, your chance of failing the test goes down. In particular,
P(Fk)=1/2*P(Fk-1), for k=2,3,4,⋯

However, the result of different exams are independent. Suppose you take the test repeatedly until you pass the test for the first time. Let X be the total number of tests you take, so Range(X)={1,2,3,⋯}.

1. Find P(X=1),P(X=2),P(X=3).
2. Find a general formula for P(X=k) for k=1,2,⋯.
3. Find the probability that you take the test more than 2 times.
4. Given that you take the test more than once, find the probability that you take the test exactly twice.

2. Relevant equations
I would say the geometric PDF: Px(k) = p(q)k-1 for integer k's, but after doing the problem I do not think it is exactly like that.
P(A|B) = P(A∩B) / P(B)

3. The attempt at a solution
This problem seemed a little too straightforward for my comfort zone, so I'm seeing if anyone will check my logic here.
P(X = 1): this is easy, it's given by the problem to be 1/2 chance that we take this test 1 time.
P(X = 2): If I solved the first part correctly, then this is equally easy - we have S2 = 1 - (1/2)(1/2) = 3/4
P(X = 3): Similarly, S3 = 7/8

This implies that P(X = k) = 1 - 1/22

Finding the probability that the test will be taken more than 2 times would be the probability that we fail on the second time which is F2 = 1/2 * F1 = 1/4

Finding the probability of 2 tests given more than 1 test is trickier for me...
Finding P(A∩B) is the intersection between X = 2 and X > 1, which is X = 2. Dividing this probability by the probability of X > 1 then is where I realize that either my initial answers were wrong or my interpretation of part 4 is wrong, but I'm leaning more towards part 1 being wrong:
P(A∩B) / P(B) = P(S2) / P(F1) = (3/4) / (1/2) = 6/4 which is impossible.

2. Sep 27, 2015

### paisiello2

The problem is a little unclear the way it is worded and you made a typo but go back to the first question:
What is the probability of P(X=2)? You gave the probability of success on the 2nd test only but we want to know the probability that you would actually end up failing the first test and then passing the 2nd test.

3. Sep 27, 2015

### whitejac

So, you're saying that rather than consider this as P(k = 2), which is what I think you're saying I did, then I should interpret P(X=2) as P(1 ∩ 2), which would be P(A)P(B)?
This is how I am understanding the problem you have with my approach - and it's one that I am feeling more certain about because the only way for X = 2 to occur would be to surpass X = 1, unlike what I said which was simply the probability of the kth term alone.

4. Sep 27, 2015

### paisiello2

Yes, I think that is what they want.

5. Sep 27, 2015

### whitejac

Doesn't this make question part 4 a little bit redundant? The only way for P(X=2) to occur would be for P(X > 1) to occur, and we are not interested in the probability of anything beyond 2.

6. Sep 27, 2015

### paisiello2

Except they added a conditional probability that you are taking it more than once.

7. Sep 27, 2015

### whitejac

That's what I thought the first time. This would be the conditional probability I mentioned in the equations section because the only requirement for (X > 1) would be the condition that F1 occurs.

8. Sep 27, 2015

### paisiello2

Yes, but you didn't get the first question right so you got the wrong answer.

9. Sep 27, 2015

### Ray Vickson

I think the best way to interpret $P(F_k)$ is as an absolute, not conditional probability. Presumably, you can write the test as many times as you want; it is just that after passing it on some trial, all further results would be ignored in assigning marks to you. It is as though you were writing additional exams just for the fun of doing it.

Of course, in the actual process you would only write Exam $k$ if you failed Exams $1 \rightarrow k-1$; but your chances of passing/failing Exam $k$ would not be influenced by previous passes or failures. So, $P(F_k) = 1/2^k$ for $k = 1,2,3, \ldots$. Naturally, we also have $P(S_k) = 1 \: - \: 1/2^k$, but that is not equal to $P(X = k)$; can you see why?

For question 3, you need to find $P(X > 2) = P(X = 3) + P(X = 4) + P(X = 5) + \cdots$. This is an infinite series. Can you see a way to reduce the problem to one having a finite summation that is readily computable?

For question 3, you want $P(X = 2|X > 1) = P(X = 2\; \cap \; X > 1) / P(X > 1) = P(X = 2)/P(X > 1)$.

10. Sep 28, 2015

### whitejac

I don't really see why yet... I read through my book more thoroughly and looked at a similar example that was the opposite:
PX(k) = 1/2k
which they proved FX(k) to be (2k - 1) / (2k) which yielded the same result as 1 - (1/2k) they did this using
x<1, FX(x) = 0
1 ≤x< 2, FX(x) = PX(1) = 1/2
2≤x< 3, FX(x) = PX(1) + PX(2) = 1/2 + 1/4 = 3/4....

This was extended to a part B that found P(2< X ≤ 5) which was
FX(5) - FX(2) = 31/32 - 3/4 = 7/32
or
PX(3) + PX(4) + PX(5) = 1/8 + 1/16 + 1/32 = 7/32

Now, if I extend this general problem and find the similarities between ours and theirs, I believe I do... But not entirely.
For my problem, then, we could define FX(x) as (1 - 1/2k), which is a nomenclature i'm still a little unsure of, but I believe means "the function for the Random Variable given x"... Then, using a definition in my book: P( a < X ≤ b) = FX(b) - FX(a):
P(X = 2) to be P(1 < x ≤ 2) = (1 - 1/22) - (1 - 1/2) ) = 1/4 => 1/2k = P(X = k) - which I suppose makes more sense numerically than my original answer because for larger value k it would seem less likely to ever happen because every other value had such an exponentially likely chance of happening as well.

11. Sep 28, 2015

### Ray Vickson

Maybe the book is doing a different example, but that does not apply in the current case. You need to think it through, not just follow the book. If X = toss number at which the first 'H' occurs, we have
$$P(\text{get heads eventually}) = P(X = 1) + P(X = 2) + P(X = 3) + \cdots$$
This is true because the first head will occur either on toss 1 or toss 2 or toss 3 .... . That means that the events $\{X = 1\}, \{X = 2 \}, \{ X = 3 \}, \ldots$ are disjoint, and so their probabilities can be safely added together to get "heads eventually'.

If you think that $P(X = 1) = 1/2, P(X = 2) = 1-(1/4), P(X = 3) = 1 - (1/8), \ldots$ what do you get for the sum? Well, already by $k = 2$ you would get a probability > 1 (because $1/2 + (1 - 1/4) = 1.25$, etc), and the sum just gets more and more above 1 as you take more terms. So, $P(X = k) = P(H_k) = 1 - 2^{-k}$ cannot possibly be right.

12. Sep 28, 2015

### whitejac

I saw why it is unable to be 1-1/2^k because, like you said, if we are to take the sum of probabilities until our numbers then we must stay below 1 and that does not. From trying to use the definition, rather than explicitely the example, I proposed that FX(2) - FX(1) would be equal to P(1 < x ≤ 2) which would give me the probability of P(x = 2). I organized this because, like you said, we have too choices repeatedly tested - pass or fail, heads or tails, and each test has a higher chance to succeed. Now, in the example and the definitions proven through the section they assigned an equation FX(x) that would be related in some way to an equation PX(x). These I thought were pass/fail designations, from the example, but if I am mistaking the nomenclature then I didn't mean to. My book never formally named F or P, they just introduced them as a method of discussing PDF's and CDF's.

Nevertheless, I can attempt to derive it logically: If we take this pass fail 'test' as P(pass eventually) = P(x=1) + P(x=2) + ... P(x=k)...
Then I am a little bit confused about this problem. It seems to me that we have 3 options on interpreting the sum of probabilities:
1) Adding up the probability of each success - we've already proven this cannot work as it grows too large.
2) Adding up the probability of each failure - this will give us a surefire way of ensuring each test occurs and it definitely sums 1. Ex: To even reach X = 9, we must fail 8 times and thus P(fail 8 times) = 1/2 + 1/4 + 1/8.... however, this does not mean X = 9, because the only way we can stop the test would be to have the test succeed.... which gives way to
3) Adding up the failures plus the final succcess - this interpretation would give me the most 'accurate' interpretation because as I see it... You add up the chance to fail until X reaches its 9th test, then you're given two choices - pass or fail. Now, failure is very small, 1/2^8, but that means that adding the chance to pass is also too large.

This final option resembles the Geometric distribution ∑pq^k, but if I put in the pass/fail constraints then I calculated the series to sum up to about 1/3 instead of 1...

13. Sep 28, 2015

### Ray Vickson

OK, you are on the right track: it is similar to a geometric distribution, but with the difference that the failure and success probabilities depend on the trial number, and so are not just constants $q = 1-p$ and $p$ at each trial. So, what do you think takes the place of $q^{k-1} p$? (Go back to square 1, and ask yourself: in the geometric case, where do the $q^k$ and the $p$ come from?)

14. Sep 28, 2015

### whitejac

Could we not say that p and q are functions? They can only be a constant probability? The geometric distribution comes from a series of Bernoulli trials, where it either passes or it does not. In my case, it passes to the tune of P(Sk) and fails to the tune of P(Fk)

15. Sep 28, 2015

### Ray Vickson

I did not ask about Bernoulli trials, etc. I just asked for an explanation of the formula q^k * p, because once you understand that, you will be able to figure out how to adapt the formula to the new situation.
So, what do you use in place of the "geometric" formula $P(X = k) = q^{k-1} p$? If you understand where the $q^{k-1}$ and $p$ come from, you ought to be able to find the new formula.

16. Sep 29, 2015

### whitejac

I must be missing something critical here... I want to view these things (the summations, p and q, the placement of k's...etc) as tools for deriving a situation mathematically but i'm not seeing the relationships or something.
p is the chance to pass, q is conversly the chance to fail.
q^k-1 is the chance to fail k-1 times and multiplied by p means that the final attempt is a success (because it has no relation to k).
Now, if we're to broaden this to allow for p relative to k then we have to account for whatever value p is k times becuse it has to be factored each time we take the test. Otherwise its chances are constant and we're left with the "parent" function. We have to assume that for every failure q, there was an opportunity p for success. This opportunity grew at a rate of 1-q, strangely enough. Charting out values of k = 1, 2... gives...

p: 1-P(Fk=1), 1-P(Fk=2) .... = 1/2, 3/4 ... 1- 1/2^k
q: 0.5, 0.5P(Fk-1).... = 1/2, 1/4 ... 1/2^k

To what means can we modify with the value of k? Simply "plugging it in" does us little as the value for p shifts every test. I proposed in my head representing p as 2^k-1/2^k, but that's equivalent to 1-q and we've already established that out of the question. Subbing in the equations to pass and fail (pq) gives us a false representation of the probability because
k = 2 => P(X=1) = 1/2 (1/2 * 1) = 1/4
k=3 => P(X=2) = 15/16 (1/2 * 1/4) = 15 / 128
...
And this gets much too small for an valid probability representation.
Bringing p^k and q^k-1 gives...
(1/4)(1/4) = 1/16
(49/64)(1/16)...
again the probability is becoming much much too small and given the fact that p is ultimately a fraction...I cannot see bringing it to the pwer of k. Futhermore, that would bring it more in line with the binomial family and that's contrary to my intentions.

17. Sep 29, 2015

### Ray Vickson

One more hint, then I quit.

Describe the sample space, in terms of sequences of F and S (fail and pass). In terms of pass/fail sequences, what does the sample point {X = 3} look like? Using that representation, how would you compute P(X = 3)? ALWAYS go back to the sample space when in doubt!

18. Sep 29, 2015

### whitejac

Okay, I'm sorry for being so obtuse at this kind of stuff. I get set in a way of looking at something and attempt to adapt it rather than step back and reevaluate it.

Considering the sample space:
We know that this is geometric in nature, and even if we didn't we know that the aim of this test is to stop at the first success. This gives us:
S = { S, FS, FFS, FFFS,...kFS}

Charting this out as tangible probabilities gives us: {(1/2), (1/2)(3/4). (1/2)(1/4)(7/8)....} = {1/2, 3/8, 7/64...}

This is the idea that if it passes then it failed all of the times before. Now, I know this to be true from the sample space "possibilities", the given probability functions, and because this is consistently approaching 1 but not quite making it. I also ran this past my professor.

The snare here is this... I'm attempting to solve this for X=K and I've found the answer but can't figure out how to write it...
The numerator is easy: (2^k - 1) => 1, 3, 7...
The denominator is a little less so...

I broke it down as follows:
2, 8, 64, 1024....
= 2^1 , 2^3, 2^6, 2^10...
= 2^1 , (2^2)(2^1) , (2^3)(2^3) , (2^4)(2^6)...
= (2^k) , (2^k)( 2^k-1) , (2^k)( 2^k-1)( 2^k-2)( 2^k-3) ..... etc

Now, I know this isn't a factorial exponent (2^k!) because it wouldn't be the same for k > 3.
Further more I know this isn't a factorial on the outside (2^k)! because then you'd get 6*5*4... and it's already wrong.

Is there a mathematical term for multiplying by the previous value rather than the subsequent integers?
My next step was to convert this into pq^k-1 format, even though it's not required for the purpose of this problem, because we discussed it to so much depth and my professor said it was possible as well as beneficial to practice the idea.

19. Sep 29, 2015

### Ray Vickson

Congratulations, you are on the right track at last----that sample space point of view really works! Anyway, if
$$Q_j = \frac{1}{2} \cdot \frac{1}{2^2} \cdots \frac{1}{2^j},$$
we have
$$P(X = k) = Q_{k-1} \left( 1 - \frac{1}{2^k} \right) = Q_{k-1} - Q_{k}.$$
This expression is handy, because from it we can get
$$\begin{array}{rcl} P(X \leq k) &=& P(X = 1) + P(X = 2) + \cdots + P(X = k)\\ &=& (1 - Q_1) + (Q_1 - Q_2) + (Q_2 - Q_3) + \cdots + (Q_{k-1} - Q_k)\\ &= &1 - Q_k \end{array}$$
because all the intermediate terms $Q_1, Q_2, \ldots, Q_{k-1}$ cancel, to leave only the first and last term.

Of course you can find $Q_j$, because
$$Q_j = \frac{1}{2^1} \frac{1}{2^2} \cdots \frac{1}{2^j} = \frac{1}{2^{1+2+ \cdots + j}},$$
and the sum of integers $1 + 2 + \cdots + j$ is well-known from any elementary algebra textbook.

Last edited: Sep 30, 2015
20. Sep 30, 2015

### whitejac

Thank you! That's an elegant way of writing pq^k.
I found P(X=k) to be (2^2k-1) / (2^k(k+1)/2).


text:
[This article was first published on Heidi's stats blog - Rbloggers, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Reproducible research is a topic that people like to talk about these days. Thinking about reproducible research and learning the important tools is what improved my work more than anything. Not in a sense that my results got better. More in a sense that my feeling about the work got better and my analyses got easier to understand for future me.

So in the following I would like to give a list of things that are helping future Heidi most:

Of course there are many other things for me that are influencing the usefullness of these tools. First of all I use R for statistical analyses and run
Linux (Ubuntu) on my computer (and servers).

This post is not meant as an instruction to how to use the tools. Introductions are linked above. I want to tell how I am connecting the tools for projects.

### How to use Literate Programming, Version Control and Makefiles together

#### knitr

As a statistician, before shipping a report, talk or paper usually many hours are spent with the cleaning and exploration of the data. In order to keep track of what I already did, I like to produce a PDF of all tables, figures and other analyses that I produce. Or at least those that are in any way usefull. That I do using knitr. With knitr I can also include my thoughts

#### Roxygen

I document (almost) all functions that I write for a project with Roxygen2. It is a fast and easy way to keep track of what your function does and what parameters are supposed to mean. And if you decide to make an R-package out of the code, it can easily be transformed into .Rd files.

#### SVN

To keep track of the changes and in order to be able to go back to old versions of my code or my latex files, I use SVN (or occasionally git). This is super nice when you work alone. It is even better when you work with other people. With SVN I control and save all my R-, Rnw-, tex- and other files that are important to keep in case my computer breaks or my office burns down. In this way it is also possible for me to work on different machines without any hassle. I just checkout the repository to the given machine.

#### Make

Imagine you write a paper or a report. You change something in the data cleaning code or in your analyses. Wouldn’t it be awesome to just do a single command and all codes that depend on your change are run and your paper/report automatically updated? It is possible by defining a simple text file (Makefile) and then just say make all. At least when you work in Linux.

#### The gold standard

My personal goal is

1. To make my work so understandable that in some years I will still understand what I did or even that I could quit a project and someone else could keep working on it without major problems.
2. To be able to go back to old thoughts and track changes that anyone did on the project.
3. To have a “flow” in my projects where I can change things somewhere in the middle of my code and with a single command (make) I can update everything that is neccesary.
4. If the worst happens and all the machines I work on brake, I want to be able to recreate everything without loosing much time.

To leave a comment for the author, please follow the link and comment on their blog: Heidi's stats blog - Rbloggers.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

# Never miss an update! Subscribe to R-bloggers to receive e-mails with the latest R posts.(You will not see this message again.)

Click here to close (This popup will not appear again)


text:
# pyspark.sql.Catalog.recoverPartitions¶

Catalog.recoverPartitions(tableName: str) → None

Recovers all the partitions of the given table and update the catalog.

Only works with a partitioned table, and not a view.


text:
# Advanced device relaxation - manual workflow¶

In this tutorial you will learn how to use a combination of “Rigid” constraints and a simple optimization routine to perform efficient and reliable relaxation of the internal coordinates in a device configuration.

Note

In the O-2018.6 release, we introduced a so-called Study Object to automate the workflows described in this tutorial. That is now the recommended way of performing device geometry optimizations in QuantumATK, and if you have access to O-2018.06 or newer, you should go to the corresponding tutorial Relaxation of devices using the OptimizeDeviceConfiguration study object. Users of previous versions can still perform the calculations using the more complex workflow, that is described here.

## Introduction¶

Relaxation of the internal coordinates of device configurations can be a tedious and time consuming task. There are two main reasons for this:

The 2-probe configuration consists of three distinct regions: Two electrodes and a central region consisting of electrode extensions and a scattering region. The geometry of each region should be separately optimized in a first-principles manner, but both electrode extensions in the central region must always match the corresponding electrode.

Step 1: Optimize the electrodes. After optimization of the electrodes, it is important to make sure that the electrode extensions in the central region match the optimized electrodes.

Step 2: Optimize the central region. The optimum length of the central region is usually not known a priori. In the case of a periodic bulk, one would minimize the stress tensor. However, the stress field is not uniquely defined for a device configuration, so another approach is needed.

In the second step, we must rely on minimizing the device total energy while explicitly preserving the internal coordinates of the electrode extensions. There are at least two different ways to do this using QuantumATK.

BRR method

After building the 2-probe configuration from fully relaxed electrodes, you extract the central region bulk and relax it with “Rigid” constraints on the atoms belonging to the electrode extensions. The device or interface is then reassembled from the relaxed bulk. This is a computationally efficient approach and will usually capture a large fraction of the required geometry relaxation. We shall refer to this method as Bulk Rigid Relaxation (BRR).

The BBR method does not calculate the total energy using a full device geometry, so the reassembled device may not be exactly in the global minimum-energy geometry. Some types of electronic structure calculations may be sensitive to this difference, others may not.

1D minimization

A somewhat more brute-force approach is to explicitly minimize the device total energy wrt. internal coordinates and the central region length, i.e., do the full 2-probe calculations and vary the central region length. This is a 1D minimization problem with relaxation of the atomic positions at each step, so we shall refer to this method as 1DMIN. We have developed a procedure for doing this as efficiently as possible with QuantumATK, but the computations are in general more demanding than with BRR. It is therefore recommended to use 1DMIN only after pre-relaxing the central region with the BRR method.

This tutorial introduces the electrode relaxation and the central region relaxation using BRR and 1DMIN methods. The methods are illustrated for relaxing a silver-gold interface using a classical potential. The classical potential is selected for computational convenience, and the same procedure could be used for another relaxation method like DFT. In general we recommend to check ATK-ForceField calculations with DFT calculations.

Tip

We have prepared a graphical explanation of the steps outlined above. You can see and download the infographic of the present tutorial at the Infographics paragraph.

## Preparations¶

First make a new project called ‘device_relaxation’. Our starting point is an unrelaxed 2-probe device configuration. It was created using the QuantumATK Interface Builder follow the tutorial at this link: Building an interface between Ag(100) and Au(111). The geometry can also be downloaded as this file: Ag_Au_interface.py.

Save the py-file in your project folder, locate the device configuration on the LabFloor, and drag-and-drop the configuration to the Builder.

The 2-probe configuration should look like the one illustrated below. It consists of a left electrode, Ag(100), the central region, and a right electrode Au(111). Use the plugin to split the device into a left electrode, right electrode and central region.

## Electrode relaxation¶

The device geometry was constructed using the interface builder (Building an interface between Ag(100) and Au(111)). Prior to building the interface the left and right crystal structures were set up using their experimental geometries. We will assume that the experimental geometry is also the lowest energy structure of the theoretical method we will use for relaxing the device geometry.

In the interface builder we choose “Strain second surface” upon selecting the surface cell. This means that the Ag(100) surface was kept in its experimental geometry, while the Au(111) surface was strained in the x,y directions upon forming the interface. When the crystal is strained in the x,y directions, it will relax in the z-direction according to the Poisson ratio of the material.

In the following we will apply this relaxation to the Au(111) crystal, by relaxing the stress in the z-direction while keeping the x, y directions fixed. The obtained relaxation of Au(111) will be transferred to the gold atoms in the central region.

Open the builder and change the name of the right electrode to Au111. Double click the geometry to make it active and send it to the Scripter. Add the following objects:

• New Calculator
• OptimizeGeometry

Now open the New Calculator object and select the ATK-ForceField “EAM_Zhou_2004” potential. Click OK.

Note

For QuantumATK-versions older than 2017, the ATK-ForceField calculator can be found under the name ATK-Classical.

Next open the Optimize Geometry object:

1. Set the force tolerance to 0.01
2. Set the stress tolerance to 0.001
3. Uncheck z in the “Constrain cell” checkboxes

Send the script to the Job Manager and run the relaxation. It should take less than a minute to execute. When finished inspect the log file.

You should now see a new object on the LabFloor called “Au111.hdf5”. It contains two configurations, corresponding to the geometry before and after the relaxation. Select the two geometries and drop them onto the Viewer. In the viewer you can hoover the mouse over the unit cell of each configuration and you will notice that the length of the C vector is 7.064 Å before the relaxation and 7.011 after the relaxation. This information could also be obtained from the log file, “Au111.log”.

The relaxation results in a small compressive strain of 0.75%. To apply this to the central region, open the Builder and double click the central region configuration.

1. Rename the stash item to “central”.
2. Select all the gold atoms by drawing a rectangle around them.
3. Open the Bulk Tools ‣ Stretch Cell plugin.
4. Enter 0.9925 along the C axis.
5. Use “Selected atoms are” stretched.
6. Click Apply.

Note

The strain is so tiny that it is not possible to see any visual change in the structure. However, you may open the “Lattice Parameters..” plugin and click the and buttons to see how the C cell vector changes when you undo and redo the stretch. Similarly, you can inspect the Coordinate List.

## Central region relaxation¶

You are now going to relax the central region only. We want the electrode extensions to be able to move freely (but rigidly) during geometry relaxation. We therefore add vacuum on both ends:

Make sure that no atoms are selected, by clicking the background of the plotting window. In Bulk Tools ‣ Stretch Cell use the C-vector slider to increase the cell length by 30% or simply type in a multiplication factor of 1.3. The cell is stretched once you click Apply.

Move to Coordinate Tools ‣ Center and uncheck A and B. Center the central region by clicking Apply.

We want to fix the atoms in the electrode extensions by including an additional layer of the central region. The additional layer will help the Device from Bulk tool in QuantumATK to identify the length of the electrode. For this purpose select the 5 outermost silver atoms. Open Selection Tools ‣ Tags. Enter Left in the box named “Click here to write a new tag...”. Press Return.

Next select the 4 outermost gold atoms, and tag them Right.

Send the configuration to the Script Generator and set up the calculation:

Add a New Calculator and an OptimizeGeometry object.

Choose the ATK-ForceField engine and the “EAM_Zhou_2004” potential for these calculations.

Make the following changes to the OptimizeGeometry object:

1. Set the force tolerance to 0.01 eV/Å.
2. Check “Save trajectory” and set the filename to central.hdf5.
3. Click “Add Constraints” to set up the Rigid constraints on the electrode extensions.
4. Check that the Left and Right tags are assigned to the correct atoms in the configuration. Do this by clicking the tags. Choose “Rigid” constraints for both tags.

Note

Choosing instead a “Fixed” constraint for one of the tagged regions while keeping the other one “Rigid”, may save BFGS relaxation steps and thereby computational time.

Send the script to the Job Manager and run the calculation. It should be very fast. Once finished, the central.hdf5 output file appears on the LabFloor.

You can use the Viewer plugin to see the trajectory of the relaxation as illustrated in the animation below.

The relaxed central region is saved in central.hdf5 as the bulk configuration with id “glD002”. This is your initial guess for the central region geometry. Drag-and-dro it on to the Builder.

Click Device Tools ‣ Device From Bulk to turn the central.hdf5 bulk into a device configuration. You can leave the settings for electrode and central region lengths at default values.

The interface configuration should now look like illustrated below. Rename the configuration to “device” by using F2 or with a right mouse-click.

## 1DMIN optimization of the interface using 2-probe calculations¶

We are now ready to search for the global minimum-energy geometry of the Ag(100)-Au(111) interface. Save the device configuration as a HDF5 file named device_in.hdf5 by choosing “Save as” after right-clicking the device Stash item.

You are now going to apply a simple SciPy minimization procedure in search for the central region length that minimizes the interface total energy. The device configuration you just saved in device_in.hdf5 is the starting point for this.

Two scripts are needed. You can download them here: device.py and optimize.py. Save both scripts in the folder where device_in.hdf5 was saved. The first script (device.py) reads in the device configuration and sets up and runs the calculation, while the second script (optimize.py) contains functions and classes used by device.py.

Note

To use the script device.py for other systems you need to edit the script and change the definition of the potential to an appropriate potential for your new system.

Both scripts should now be available in the QuantumATK “Project Files” list. Drag-and-drop device.py onto Jobs. The Job Manager will open. Click to run the script.

Once finished (it’s fast!), a window with the figure below will pop up. The blue points are total energies of relaxed interface configurations. The initial cell length is 28.51 Å and the other blue points are trial guesses made by the algorithm. The red point indicates the optimum central region length which is 28.58 Å. The optimization routine simply maps out the illustrated potential-energy surface and finds its minimum!

The job will finish when you close the pop-up window.

The relaxed configurations are saved in the output file device_out.hdf5 along with the total energies. The configuration corresponding to the red point (the final device or interface geometry) is the latest device configuration (the one with largest id number).

Using the Viewer plugin, the minimum-energy interface configuration looks like this:

In this case the change to the structure after the 1DMIN relaxation is minute compared to the BRR relaxation. Often the BRR will be sufficiently accurate and there is no need for the 1DMIN step. The 1DMIN approach is only needed for extreme accuracy.

## Infographics¶

Download Infographics in pdf version. Click the button below to see the infographics.