text:
# Why does instant recoil after a punch deliver a stronger blow?

## Recommended Posts

I practice kenpo karate, and in our technqiues for punching, we are told to deliver a punch but then as soon as it hits the target, recoil it back to yourself as fast as possible.

I've tested this out numerous times on a punching bag (hi-tech right?), and every time I throw a punch and then recoil it, I move the punching bag further than just throwing a

one direction punch.

-Just curious

~EE

##### Share on other sites

I practice kenpo karate, and in our technqiues for punching, we are told to deliver a punch but then as soon as it hits the target, recoil it back to yourself as fast as possible.

I've tested this out numerous times on a punching bag (hi-tech right?), and every time I throw a punch and then recoil it, I move the punching bag further than just throwing a

one direction punch.

-Just curious

~EE

That could the fault of the device that is measuring the punch. (That was not supposed to be a joke but is now!)

Edited by Robittybob1

##### Share on other sites

That could the fault of the device that is measuring the punch. (That was not supposed to be a joke but is now!)

That hurt more than getting punched in the face.

##### Share on other sites

That hurt more than getting punched in the face.

Now you are being a softy.

We were taught the same, to jab real fast, I thought it was to have your hands back in the defensive position as quick as possible.

##### Share on other sites

Now you are being a softy.

We were taught the same, to jab real fast, I thought it was to have your hands back in the defensive position as quick as possible.

Well ya thats the point, but it seems to be a stronger punch as well. I was wondering if it was something similar to how a bullet is spun down a barrel, making it more accurate

##### Share on other sites

It may have something to do with the speed of the punch. Without actually seeing you punch I would say that you are hitting faster when you recoil than when you throw your 'one direction punch'. Punching power is more appropriately mesured by how much you compress the bag (inelastic collision) than it is how far you move it (elastic collision). The more the bag moves the less impact the bag feels as energy that is meant for the bag is lost in the bags motion.

##### Share on other sites

Well ya thats the point, but it seems to be a stronger punch as well. I was wondering if it was something similar to how a bullet is spun down a barrel, making it more accurate

Rifling leverages conservation of angular momentum, which does not appear to be an issue with a jab.

———

If this were a scattering interaction, the difference would be between a completely inelastic collision and an elastic one (or close to that). In an elastic collision, the final speed of the target is twice that as in a completely inelastic collision (the recoil means a greater transfer of momentum) (v2f = 2m1/(m1+m2) vs m1/(m1+m2) )

The problem is that this probably can't be modeled that way, since you are actively pulling the fist back, rather than it happening as a result of the dynamics of the collision. But if it's true that the collision is closer to elastic in one case, then you will transfer more energy.

edit: as I think about this more, I think recoil is the wrong metric. Inelastic collisions transfer less kinetic energy to the target, but that lost KE appears as heat, sound and deformation, much of which would count as damage of a blow. There's more to it than recoil.

##### Share on other sites

The way this was explained to me (read: I'm not certain this is is accurate) is that punching and leaving your fist there allows more of the energy to come back and be reabsorbed by your arm. Punching and withdrawing your fist rapidly ensures the energy stays with the target in full.

Whether or not that's true, the act of snapping your punch or tensing your muscles in such a way to allow you to pull back your fist from the target probably allows just that little bit of extra "oomph" to be added to the strike... Just those few extra muscles to add to the overall impact.

##### Share on other sites

The way this was explained to me (read: I'm not certain this is is accurate) is that punching and leaving your fist there allows more of the energy to come back and be reabsorbed by your arm. Punching and withdrawing your fist rapidly ensures the energy stays with the target in full.

Whether or not that's true, the act of snapping your punch or tensing your muscles in such a way to allow you to pull back your fist from the target probably allows just that little bit of extra "oomph" to be added to the strike... Just those few extra muscles to add to the overall impact.

There seems to be a high point as you pull back that your fist/arm/body is at its most rigid, somehow allowing energy to conduct better.

##### Share on other sites

I practice kenpo karate, and in our technqiues for punching, we are told to deliver a punch but then as soon as it hits the target, recoil it back to yourself as fast as possible.

I've tested this out numerous times on a punching bag (hi-tech right?), and every time I throw a punch and then recoil it, I move the punching bag further than just throwing a one direction punch.

-Just curious

~EE

I suggest taking some actual measurements rather than trusting your impressions. Put a scale with vertical markings on one side of the bag and a video camera opposite and on the other side of the bag and then deliver your blows from 90º relative to the camera's line-of-sight. Make multiple recordings and carefully review the results and report back.

There is no physical reason to expect a quickly withdrawn blow delivers more energy, though there is a physical reason that a quicker blow does.

The amount of energy delivered to the target by the hammer-blow is equivalent to one half the mass of the head times the square of the head's speed at the time of impact . While the energy delivered to the target increases linearly with mass, it increases quadratically with the speed (see the effect of the handle, below). ...

##### Share on other sites

Why not make a real experiment from this?

Get Force Gauge Analog/Digital Device.
http://www.industrial-needs.com/technical-data/force-gauge-gram-division.htm

Connect wire to it, then to bottom of bag.

Somebody with digital camera will have to record device while hitting it.

##### Share on other sites

Why not make a real experiment from this?

Get Force Gauge Analog/Digital Device.

Connect wire to it, then to bottom of bag.

Somebody with digital camera will have to record device while hitting it.

Are you suggesting my suggestion was not a real experiment!

That's a good addition however. I would also suggest having the camera positioned so that the entire puncher is in the frame so that it's possible to gauge the speed of the punch and any differences in stance or arm position.

##### Share on other sites

Are you suggesting my suggestion was not a real experiment!

Not at all my friend..

But you were suggesting using camera to record movement of bag (like in pendulum), isn't?

The more powerful hit, bigger deviation from natural vertical direction of bag.

It would be hard to get good precision. And better to have 1000 FPS camera with it.

BTW, I found that GoPro Hero 4 for just $500 can record at 240 FPS: Not bad, and has also quite nice configurable interval recording mode/time lapse 1s-60s. 120 FPS @ Full HD 30 FPS @ 4K UHD That's a good addition however. I would also suggest having the camera positioned so that the entire puncher is in the frame so that it's possible to gauge the speed of the punch and any differences in stance or arm position. Then we need two cameras. From my own experience of recording of experiments (and I am using Full HD) we need close look up to to see digits, not to mention analog arrow.. Edited by Sensei #### Share this post ##### Link to post ##### Share on other sites Not at all my friend.. But you were suggesting using camera to record movement of bag (like in pendulum), isn't? The more powerful hit, bigger deviation from natural vertical direction of bag. It would be hard to get good precision. And better to have 1000 FPS camera with it. Just joshin' ya Sensei. But yes, the more energy directed into the bag the further it swings out. BTW, I found that GoPro Hero 4 for just$500 can record at 240 FPS:

Not bad, and has also quite nice configurable interval recording mode/time lapse 1s-60s.

120 FPS @ Full HD

30 FPS @ 4K UHD

You find the nicest equipment! A bit off-topic, but if I didn't say elsewhere, I ended up getting a Coolpix 830L. It doesn't have the interval record but it has great macro, panaorama vertical & horizontal, and HD video. Thnx for putting me onto it.

Then we need two cameras. From my own experience of recording devices on digital camera (and I am using Full HD) we need close look up to to see digits, not to mention analog arrow..

Probably add at least a 3rd camera directly behind the bag looking toward the puncher to detect any angle in the swing other than directly away from the puncher.

##### Share on other sites

You find the nicest equipment! A bit off-topic, but if I didn't say elsewhere, I ended up getting a Coolpix 830L. It doesn't have the interval record but it has great macro, panaorama vertical & horizontal, and HD video. Thnx for putting me onto it.

Do you have wifi in your model?

If yes, then you could try what I can do in S6600, and control it by smartphone:

Then write code in C/C++/Java(?) in smartphone app to take photo on demand in schedule that you want!

Edited by Sensei

##### Share on other sites

Do you have wifi in your model?

If yes, then you could try what I can do in S6600, and control it by smartphone:

Then write code in C/C++/Java(?) in smartphone app to take photo on demand in schedule that you want!

No WIFI. It's L830; not 830L as I said. Here's the specs: >COOLPIXL830

Back on topic, besides an angled swing showing up as energy not straight onto the target, heavy bags are usually suspended by chains and any upward force in a strike will make the bag jump up and so swing out less.

All-in-all, I think the "instant recoil after a punch delivers a stronger blow" idea is a myth. Happy punching!

##### Share on other sites

...

All-in-all, I think the "instant recoil after a punch delivers a stronger blow" idea is a myth. Happy punching!

My cynical intuition agrees with you. There seems to be over-reliance on received wisdom and dogma in much of the eastern martial arts.

One of my major objections would be the taught technique of most striking sports where regain of posture and defence is not paramount. In any form of martial art the greatest part is not the power of the strike but the ability to be in control and command of both yourself and your opponent - you do not do this by over-extending, nor by being overbalanced by losing control of a strike. The ability to strike with maximum force is sacrificed for composure and to maintain stance that allows defence/counter/renewal of offence; the snapped back punch is just one method of ensuring that your attack is not a be-all-end-all gamble. However in many non-martial sports the object is to impart maximum controlled force in a one off situation; teeing off on the long par four, the straight six in cricket, taking three points from a penalty inside your own half in Rugby - and all of these rely on the long, controlled and fluid follow through. In these sports we are taught that to dampen the shot, we relax the wrist, collapse the elbow, pull the follow through; but that for a power shot the follow-through is crucial.

And a lovely excuse to post a video of the amazing Chris Gayle in striking form

##### Share on other sites

My cynical intuition agrees with you. There seems to be over-reliance on received wisdom and dogma in much of the eastern martial arts.

One of my major objections would be the taught technique of most striking sports where regain of posture and defence is not paramount. In any form of martial art the greatest part is not the power of the strike but the ability to be in control and command of both yourself and your opponent - you do not do this by over-extending, nor by being overbalanced by losing control of a strike. The ability to strike with maximum force is sacrificed for composure and to maintain stance that allows defence/counter/renewal of offence; the snapped back punch is just one method of ensuring that your attack is not a be-all-end-all gamble. However in many non-martial sports the object is to impart maximum controlled force in a one off situation; teeing off on the long par four, the straight six in cricket, taking three points from a penalty inside your own half in Rugby - and all of these rely on the long, controlled and fluid follow through. In these sports we are taught that to dampen the shot, we relax the wrist, collapse the elbow, pull the follow through; but that for a power shot the follow-through is crucial.

The difference with a Karate strike is that the target is not propelled with just a momentary contact unless the strike is recoiled. then, you have a similar scenario.

##### Share on other sites

I don't know if they still sell them but you can actually buy pads and punchinng bags that have force sensors built into them if you can afford them.

##### Share on other sites

The difference with a Karate strike is that the target is not propelled with just a momentary contact unless the strike is recoiled. then, you have a similar scenario.

I think - as this is in physics rather than the lounge - to be more precise in what we say and what we are seeking. Contacts are never instantaneous (and I note that you didn't say they were) - they last whilst compression waves travel back and forth and look amazing in slo-mo. It might be worth checking if the slow-mo guys have videos that could help in judging the action.

What does a strong punch do? What makes a karate blow powerful? Let us try to separate physics and mechanics from mysticism and the sporting/point-scoring ideas of much of modern martial arts

##### Share on other sites

When you guys have all finished discussing force gauges cameras and punch bags, perhaps you could address the difficult part.

How do you propose to eliminate the psychological factor of hitting the bag harder because you think you can hit it harder?

##### Share on other sites

When you guys have all finished discussing force gauges cameras and punch bags, perhaps you could address the difficult part.

How do you propose to eliminate the psychological factor of hitting the bag harder because you think you can hit it harder?

I bow to your superior intellect in explicating the matter for us.

##### Share on other sites

I bow to your superior intellect in explicating the matter for us.

Which is the implication in #2. It is the brain controlling the arm movement, so when you decide to jab quickly you are less concerned about reaching. So the arm moves differently in both cases.

##### Share on other sites

The way this was explained to me (read: I'm not certain this is is accurate) is that punching and leaving your fist there allows more of the energy to come back and be reabsorbed by your arm. Punching and withdrawing your fist rapidly ensures the energy stays with the target in full.

This is what I was looking for, relative to an explanation. Now to find out if this actually holds any truth.

##### Share on other sites

The way this was explained to me (read: I'm not certain this is accurate) is that punching and leaving your fist there allows more of the energy to come back and be reabsorbed by your arm. Punching and withdrawing your fist rapidly ensures the energy stays with the target in full.

This is what I was looking for, relative to an explanation. Now to find out if this actually holds any truth.

Since the energy delivered to the target is a function of the mass and speed of/in the punch, then that explanation strikes me as clearly inaccurate. Energy 'reabsorbed by your arm' makes no physical sense here. The puncher does not gain any mass which leaves speed as a factor and how does one absorb speed? Edited by Acme

## Create an account

Register a new account


text:
# Conditional expectation $E[\min{\{X,Y\}}|Y^{2}]$

Let $X$ and $Y$ be an independent random variables with exponential distribution $Exp(1)$.

Compute $E[\min{\{X,Y\}}|Y^{2}]$.

My attempt:

I. If we want to compute $E[\min{\{X,Y\}}|Y]$ may be sufficient to use the following formula:

$$E[f(X,Y)|Y]=\int\limits_{-\infty}^{+\infty}f(x,Y)f_{X}(x)dx$$

Is this correct? Is there any other simpler method to compute this conditional expectation?

II. What about $E[\min{\{X,Y\}}|Y^{2}]$? First thing which comes to my mind is that we should find distribution of vector $(\min{\{X,Y\}},Y^{2})$ and $Y^{2}$ and then use appropriate formula. Am I right? How to find distribution of vector $(\min{\{X,Y\}},Y^{2})$? Is there any simpler method to compute this conditional expecation?

• HINT: if $X \sim Exp(\lambda_x)$ and $Y \sim Exp(\lambda_y)$, then $Z = \min\{X,Y\} \sim Exp(\lambda_x + \lambda_y)$ – Guilherme Thompson Dec 14 '15 at 9:11
• @GuilhermeThompson ?? "HINT" is not useful and "HINT2" is wrong. – Did Mar 19 '16 at 10:43

1. Yes, you are correct. So we have

$$\Bbb{E}[\min\{X, Y\} | Y] = \int_{0}^{\infty} \min\{x, Y\} \, e^{-x} \, dx = 1 - e^{-Y}.$$

2. Notice that both $Y$ and $Y^2$ generate the same $\sigma$-algebra: $\sigma(Y) = \sigma(Y^2)$. So we have

$$\Bbb{E}[f(X, Y)|Y^2] = \Bbb{E}[f(X, Y)|\sigma(Y^2)] = \Bbb{E}[f(X, Y)|\sigma(Y)] = \Bbb{E}[f(X, Y)|Y].$$

If you have not learned the axiomatic approach to the probability theory, ignore this explanation. We have an alternative, more computational explanation as follows: Let $Z = Y^2$. Then

\begin{align*} \Bbb{E}[f(X, Y) | Y^2] &= \Bbb{E}[f(X, \sqrt{Z}) | Z] \\ &= \int_{0}^{\infty} f(x, \sqrt{Z}) \, e^{-x} \, dx \\ &= \int_{0}^{\infty} f(x, Y) \, e^{-x} \, dx \\ &= \Bbb{E}[f(X, Y) | Y]. \end{align*}

• @Leon, Yes, any exponential r.v. is non-negative with probability 1. – Sangchul Lee Dec 14 '15 at 9:50
• I am stuck at a similar problem, would you mind adding a comment on what you did to compute the integral? i.e. $$\int_0^\infty \min \lbrace x,Y \rbrace e^{-x} dx = ?$$ – Spaced Oct 1 '16 at 18:34
• @Spaced, Have you tried splitting the integral into two parts according to whether $x < Y$ or $x > Y$? – Sangchul Lee Oct 1 '16 at 19:58
• Thanks a lot that did the trick! (the best thing being that it's not even a trick, more like an application of the indicator function) – Spaced Oct 1 '16 at 20:49

format: plain


text:
## Winter 2016 Final Exam Q2A [ENDORSED]

$\Delta G^{\circ}= \Delta H^{\circ} - T \Delta S^{\circ}$

$\Delta G^{\circ}= -RT\ln K$

$\Delta G^{\circ}= \sum \Delta G_{f}^{\circ}(products) - \sum \Delta G_{f}^{\circ}(reactants)$

Absolute Zero 1O
Posts: 14
Joined: Wed Sep 21, 2016 2:56 pm

### Winter 2016 Final Exam Q2A

Just a quick question, are partial pressures always supposed to be in bars? Once we have K=2x/x^2, how did we know that the unit for x would be bar, and not atm for example?

Claudette_Contr_3I
Posts: 23
Joined: Wed Sep 21, 2016 2:57 pm

### Re: Winter 2016 Final Exam Q2A

I think this is because bar may be the standard unit for partial pressures. Whether you need the partial pressure in bar or in atm form relies on the question.

Stephanie Demo 2N
Posts: 11
Joined: Thu Sep 29, 2016 3:00 am

### Re: Winter 2016 Final Exam Q2A

Last quarter in 14A, Professor Li said that it is common practice to use bar as the unit of pressure, but I'm not sure if Dr. Lavelle specified using bar or atm. Just be sure that you are consistent with units when doing calculations so that you don't lose points for conversions.

Aliya Habib 1L
Posts: 20
Joined: Fri Sep 25, 2015 3:00 am

### Re: Winter 2016 Final Exam Q2A

I think you just work with the units you are given in the question.

Raul Hernandez
Posts: 28
Joined: Fri Jun 17, 2016 11:28 am

### Re: Winter 2016 Final Exam Q2A

Kind of unrelated question but: why is it 2X on the top? I didn't quite get that in the review session. Do you mind explaining that part?

Jennifer Ngo 1M
Posts: 23
Joined: Sat Jul 09, 2016 3:00 am

### Re: Winter 2016 Final Exam Q2A

Raul Hernandez wrote:Kind of unrelated question but: why is it 2X on the top? I didn't quite get that in the review session. Do you mind explaining that part?

The reactants are half of the products aka the products are double the reactants. He made the concentration of the reactants "x," so products would be "2x." You also could have made the products "x" and reactants "x/2"

Chem_Mod
Posts: 17997
Joined: Thu Aug 04, 2011 1:53 pm
Has upvoted: 418 times

### Re: Winter 2016 Final Exam Q2A  [ENDORSED]

Bar is the SI unit for pressure (and very close to ATM).

The question stated, "Calculate the partial pressures in bar."

1.00 atm = 1.01 bar Very similar, which is why we often use atm even though bar is the official SI unit for pressure (as I mentioned several times in Chem 14A and in my 2 hour review today.

This was discussed in detail in my review.

Kiara1F
Posts: 64
Joined: Fri Sep 25, 2015 3:00 am

### Re: Winter 2016 Final Exam Q2A

why is there an x^2 on the bottom?
Is it because there are 2 reactants?
It just threw me off because I was thinking in terms of only squaring the number if they had coefficients in front of them but they don't. But there are 2 sets of reactants....


text:
# Inverting the infinite matrix $+\mathbf{I}$ with entries $\mathbf{P}_{ij}={i-1\choose j-1}$ [closed]

Let $\mathbf{P}$ denote the "infinite matrix"

$$\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & \dots \\ 1 & 1 & 0 & 0 & \dots \\ 1 & 2 & 1 & 0 & \dots \\ 1 & 3 & 3 & 1 & \dots \\ \vdots & \vdots & \vdots & \vdots & \ddots \\ \end{array} \right]$$

with entries $\mathbf{P}_{ij} = \dbinom{i-1}{j-1}$ and let $\mathbf{I}$ denote the "infinite identity matrix." Compute the inverse of $\mathbf{P} + \mathbf{I}$.

This was not the initial attempt. I couldn't think of anything at first. But after some nudges, I tried to compute the $n\times n$ matrix $\mathbf{P}_n$. Now my first observation was $\det \mathbf{P}_n=1$. Now if I could show that the invertibility of $\mathbf{P}_n$ would be efficient. So we expand $\mathbf{P}_{n+1}$ by the last row, then it is obvious that $\det \mathbf{P}_{n+1}=\det \mathbf{P}_n=1$. So invertibility is meaningful. But when I inverted for small values, I couldn't find any pattern. I can't think of a method to cook up the solution. Can someone help me? I see my method of thinking should have been presented and I apologise. I will add them later on.

How to compute an inverse of an infinite matrix? And even if I can, what to do with it? Thanks for any help.

## closed as off-topic by heropup, colormegone, Najib Idrissi, Claude Leibovici, Jean-Claude ArbautAug 24 '14 at 9:30

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question is missing context or other details: Please improve the question by providing additional context, which ideally includes your thoughts on the problem and any attempts you have made to solve it. This information helps others identify where you have difficulties and helps them write answers appropriate to your experience level." – heropup, colormegone, Najib Idrissi, Claude Leibovici, Jean-Claude Arbaut
If this question can be reworded to fit the rules in the help center, please edit the question.

• Why do you want to compute the inverse? Who told you there is something you can do with it? – Jonas Meyer Aug 21 '14 at 15:12
• My teacher gave me this problem. He said he has done it. I hope this satisfies your queries? @JonasMeyer Thanks. Now could you help me out here? Thanks again. – shadow10 Aug 21 '14 at 15:16
• It's possible the teacher just wants you to discover some pattern in the inverse, then prove that this pattern holds for the whole matrix. It looks like this matrix contains Pascal's triangle, although I don't know for sure. – rschwieb Aug 21 '14 at 15:16
• It looks like $P_{ij}^{-1}=(-1)^{i+j}P_{ij}$. – g.kov Aug 21 '14 at 19:32
• @shadow10 It seems like there was a rather obvious pattern emerging right at the start, which I think is what g.kov was mentioning. Did it break down on some row? As I mentioned, there's no point in talking about determinants for an infinite matrix. If you find a formula for the entries of the inverse, then you can just check the matrix multiplication directly to see if it is something like the Kronecker delta function. – rschwieb Aug 22 '14 at 12:29

## 2 Answers

For now, let's focus on inverting $\mathbf{P}$. (The method below should work for $\mathbf{I}+\mathbf{P}$ as well but I didn't want to start there). First, note that the elements of $\mathbf{P}^{-1}$ satisfy $$(\mathbf{P}\cdot\mathbf{P}^{-1})_{ik}=\sum_{j} \mathbf{P}_{ij}(\mathbf{P}^{-1})_{jk} =\sum_j \binom{i-1}{j-1}(\mathbf{P}^{-1})_{jk}=\delta_{ik}.$$ Note that this is still a linear algebra problem, albeit with an infinite number of variables and constraints.

I'll attack it with a generating function approach: Multiplying the LHS by $x^i$ and summing over all integers produces

\begin{align} \sum_{ij} \binom{i-1}{j-1}(\mathbf{P}^{-1})_{jk}x^i=\sum_j (\mathbf{P}^{-1})_{jk}\sum_i\binom{i-1}{j-1}x^i = \sum_{j=1}^\infty (\mathbf{P}^{-1})_{jk}\left(\frac{x}{1-x}\right)^j. \end{align} To justify the last equality, shift the index of summation of $i\mapsto i+j$: $$\sum_i\binom{i+j-1}{j-1}x^{i+j}=x^j\cdot \sum_i\binom{i+j-1}{i}x^{i}=\dfrac{x^j}{(1-x)^j}$$ since the last summation is a (negative) binomial series. If we repeat this on the RHS we simply get $x^k$ since the Kronecker delta kills the rest of the terms. We then let $y=\dfrac{x}{1-x}$ and equate the RHS and LHS to obtain

$$\sum_{j=1}^\infty (\mathbf{P}^{-1})_{jk}y^k =\left(\frac{y}{1+y}\right)^k=(-1)^k\left[\frac{(-y)}{1-(-y)}\right]^k=(-1)^k\cdot \sum_j\binom{j-1}{k-1}(-y)^j$$ with the last equality following from the prior equation. Identifying coefficients on both sides then finally gives $\boxed{(\mathbf{P}^{-1})_{jk}=(-1)^{j+k} \binom{j-1}{k-1}}$. Comparing with our original equation, this implies the summation $\sum_j (-1)^{j+k} \binom{i-1}{j-1}\binom{j-1}{k-1}=\delta_{ik}$; this almost certainly admits a counting proof via inclusion-exclusion.

So we may take $\mathbf{P}^{-1}$ as known, and can now focus on $(\mathbf{I}+\mathbf{P})^{-1}$. Something like the binomial inverse theorem should come in handy; I'll see if I can find a simple route.

Let $V$ be the vector space of polynomials (or polynomial functions if you prefer that). Let us look at the transpose of $P$ instead of $P$. We see immediately that $P^T$ is the matrix of the linear transformation $S:p(x)\mapsto p(x+1)\in GL(V)$ with respect to the natural basis $\{1,x,x^2,\ldots\}$. Therefore the inverse of $P^T$ is the transpose of the matrix of the inverse transformation $S^{-1}:p(x)\mapsto p(x-1)$. This proves the observation/conjecture in a comment by g.kov: $P^{-1}_{ij}=(-1)^{i+j}P_{ij}$.

The matrix $I+P^T$ corresponds to the transformation \begin{aligned} R:p(x)\mapsto & p(x)+p(x+1)\\ =&p(x)+p(x)+Dp(x)+\frac{D^2}2p(x)+\cdots\\ =&(1+e^D)p(x). \end{aligned} Therefore its inverse is $$R^{-1}=\frac1{1+e^D}=\frac12+\frac12\cdot\frac{1-e^D}{1+e^D}=\frac12(1+\tanh\frac D2).$$ The Taylor series of the hyperbolic tangent is too scary for me, so may be I stop.

• You might find the two links I posted above in comments earlier interesting for comparison with that last result. (I'll admit, once I saw those two papers I gave up on there being a simple approach.) – Semiclassical Aug 24 '14 at 23:41
• @Semiclassical: Thanks. I did miss those links. My main point was to give a simple argument for $P^{-1}$. – Jyrki Lahtonen Aug 25 '14 at 4:24
• The series begins with $1/2 + 1/4 D - 1/48 D^3 + 1/480 D^5 - 17/80640 D^7 + 31/1451520 D^9$ $- 691/319334400 D^{11} +...$. The numbers $17,31,691$ in connection with the index in steps of 2 should ring some bell in direction of Bernoulli-numbers, or even better, of Euler numbers (see wikipedia). The remaining seems to be not too difficult... – Gottfried Helms Sep 27 '14 at 1:24


text:
Evolutionary Dynamics¶

We will now consider how Game Theory can be used to study evolutionary processes. The main difference is that we now consider not two player games but game with an infinite population. The strategies will make up a dynamic population that changes over time.

Reproduction¶

Consider a simple model of populatio growth: let $x(t)$ denote the size of the population at time $t$ and let as assume that the rate of growth is $a$ per population size:

$$\frac{dx}{dt}=ax$$

Note that from here on we will refer to this rate as a fitness.

The solution of this differential equation is:

$$x(t)=x_0e^{at}\text{ where }x_0=x(0)$$
In [10]:
import sympy as sym
sym.init_printing()
x = sym.Function('x')
t, a = sym.symbols('t, a')
sym.dsolve(sym.Derivative(x(t), t) - a * x(t), x(t))

Out[10]:
$$x{\left (t \right )} = C_{1} e^{a t}$$

(This is exponential growth.)

We can also use scipy to solve this differential equation numerically (relevant for more complex dynamics):

In [11]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from scipy.integrate import odeint

t = np.linspace(0, 10, 100)  # Obtain 100 time points

def dx(x, t, a):
"""Define the derivate of x"""
return a * x

If $a=10>0$:

In [12]:
a = 10
xs = odeint(func=dx, y0=1, t=t, args=(a,))
plt.plot(xs);

If $a=-10<0$:

In [13]:
a = -10
xs = odeint(func=dx, y0=1, t=t, args=(a,))
plt.plot(xs);

Selection¶

Reproduction alone is not enough to study evolutionary processes. Let us consider a population of made up of two types of individuals:

• $x(t)$ denotes the first type;
• $y(t)$ denotes the second type.

Let us assume the same expressions for the as before:

$$\frac{dx}{dt}=ax\qquad\frac{dy}{dt}=bx$$

both these population will increase or decrease independantly so there's not much of interest there but if we introduce the following:

$$\rho(t) = \frac{x(t)}{y(t)}$$

then $\lim_{t\to\infty}\rho(t)$ indicates which type takes over the population over time.

We have:

$$\frac{d\rho}{dt} = \frac{\frac{dx}{dt}y - \frac{dy}{dt}x}{y ^ 2} = \frac{xy(a - b)}{y^2}$$

which gives:

$$\frac{d\rho}{dt} = (a-b)\rho$$

which has solution (this is just the same differential equation as the previous section):

$$\rho(t) = \rho_0e^{(a-b)t}\text{ where }\rho_0=\rho(0)$$

note that even if both population grow, but one grows faster than the other (eg $a > b$) then the overall population will grow but one will take over:

In [14]:
def drho(rho, t, a, b):
"""Define the derivate of x"""
return (a - b) * rho

a, b = 10, 5
rhos = odeint(func=drho, y0=1, t=t, args=(a, b))
plt.plot(rhos);

Selection with constant population size¶

Let us consider the case of $x(t) + y(t)=1$: so the case of a constant population size (choosing a constant of 1 is just a question of scale). For this to be possible, the rates need to reduced:

$$\frac{dx}{dt}=x(a - \phi)\qquad\frac{dy}{dt}=y(b - \phi)$$

because $x(t) + y(t)=1$:

$$\frac{dx}{dt} + \frac{dy}{dt} = 0$$

also:

$$\frac{dx}{dt} + \frac{dy}{dt} = ax + by - \phi(x + y)= ax + by - \phi$$

thus $\phi=ax+by$ (this corresponds to the average of the fitness).

Substituting $y=1-x$ we have:

$$\frac{dx}{dt}=x(a - ax-b(1-x))=x(a(1 - x)-b(1-x))$$

giving:

$$\frac{dx}{dt}=x(a-b)(1-x)$$

We do not need to solve this differential equation. There are two stable points:

• $x=0$: no population of first type: no change
• $x=1$: no population of second type: no change

Also:

• $a=b$: if both types have the same fitness: no change
In [15]:
def dxy(xy, t, a, b):
"""
Define the derivate of x and y.
It takes xy as a vector
"""
x, y = xy
phi = a * x + b * y
return x * (a - phi), y * (b - phi)

a, b = 10, 5
xys = odeint(func=dxy, y0=[.5, .5], t=t, args=(a, b))
plt.plot(xys);

In [16]:
a, b = 10, 5
xys = odeint(func=dxy, y0=[1, 0], t=t, args=(a, b))
plt.plot(xys);

In [17]:
a, b = 10, 5
xys = odeint(func=dxy, y0=[0, 1], t=t, args=(a, b))
plt.plot(xys);

In [18]:
a, b = 5, 5
xys = odeint(func=dxy, y0=[.5, .5], t=t, args=(a, b))
plt.plot(xys);


text:
## Algebra 2 (1st Edition)

Published by McDougal Littell

# Chapter 12 Sequences and Series - 12.3 Analyze Geometric Sequences and Series - 12.3 Exercises - Problem Solving - Page 816: 60b

#### Answer

$b_n= (\dfrac{8}{9})^n S\\b_{12}= (\dfrac{8}{9})^{12} S$

#### Work Step by Step

Let $S$ be the area of the initial square. Here, we have $b_1= S-\dfrac{S}{9}$ and $b_2= \dfrac{8S}{9}-(8) \dfrac{S}{(9)^2}=\dfrac{64S}{81}$ Likewise, we have $b_n= (\dfrac{8}{9})^n S$ Hence, $b_{12}= (\dfrac{8}{9})^{12} S$

After you claim an answer you’ll have 24 hours to send in a draft. An editor will review the submission and either publish your submission or provide feedback.


text:
favorite We need a little bit of your help to keep things running, click on this banner to learn more
Competitions

# Number Theory

For the given positive integer n find the number of integers m, such that 1mn, gcd(m, n) ≠ 1 and gcd(m, n) ≠ m. gcd is an abbreviation for "greatest common divisor".

#### Input

Each line contains one positive integer n (0 < n < 231).

#### Output

For each number n print the number of required integers m.

Time limit 1 second
Memory limit 128 MiB
Input example #1
1
6
10
2147000000

Output example #1
0
1
3
1340599805


text:
# 1.6. Processing Elmer/Ice output

## 1.6.1. Felicity Holmes

### 1.6.1.1. General Project Description

For my PhD project, I use an ice-sheet model called Elmer/Ice. Specifcally, I am simulating the glacier Kronebreen, Svalbard for the time period 2016 - 2017 and comparing the output to observational data in order to evaluate the applicability of Elmer/Ice in this location. The long-term goal is to use the tuned model to simulate the development of the glacier until 2100.

After each simulation, I have the following output:

1. A .csv file with point data for various variables (e.g. velocity, pressure, elevation) at each timestep

2. A .dat file for each partition containing the mean velocity value after each timestep. Of particular interest is partition 0, which covers the front (and most dynamically exciting) part of the glacier.

For the project, I want to explore this data and build up a processing chain that I can use to analyse the output. My current simulations in Elmer/Ice are not final as I am still tuning parameters etc, but I use my ‘best yet’ output for this project. This means that any conclusions drawn from the analysis are only preliminary, but provide a framework for exploration and analysis of the output.

Specifically, my aims are: 1) To convert the point data into a raster than I can then compare to observations. This would allow me to identify whether the model is performing well and highlight any areas of particulary good or bad performance.

1. To look at the time series of mean velocities and understand how important surface mass balance (and its constituent components) are with regards to determining velocities. In order to do this, I want to explore the utility of a few techniques such as dynamic time warping and a random forest regressor.

I will split up this report into parts 1 and 2; part 1 will deal with the point data (mainly geocomputation) and part 2 will deal with the data from the frontal partition (geocomputation and modelling).

### 1.6.1.2. Part 1: Point data exploration

Although it is possible to export a .csv file for every timestep, I have opted to focus on the final timestep for this project. The initial conditions for the glacier are derived from observations, so it is of interest to see if the modelled and observed development of the glacier are similar. Although it would be of interest to conduct this analysis at regular intervals (e.g. once a month), the output that I am using for this project is preliminary and only spans 52 days. Therefore, this project will only compare the final model output to the appropriate observational data but forms a basis for more comprehensive analysis in the future.

The .csv data is 3D with data being provided for the entire depth of the glacier. However, I am primarily interested in the surface values as the observational data only shows surface properties. The first steps are therefore to:

1. Extract all the data when ‘depth < 5’. I am using depth < 5 rather than depth == 0 as the surface nodes do not always lie precisely at 0m depth due to mesh evolution

2. Investigate the available variables and the number of data points

[13]:

%%bash

# Extract surface data. Using depth < 5 rather than depth == 0 as some surface nodes
# are not exactly 0
awk '{if ($0<5) print}' "./inputs/calving_elmer_output.txt" > ./outputs/depth0_calvingVariables.txt head -n 1 ./outputs/depth0_calvingVariables.txt wc -l ./outputs/depth0_calvingVariables.txt tail -n 1 ./outputs/depth0_calvingVariables.txt | grep -o " " | wc -l  "depth" "temperature homologous" "smb" "velocity:0" "velocity:1" "velocity:2" "frontadvance:0" "frontadvance:1" "frontadvance:2" "calving:0" "calving:1" "calving:2" "Points:0" "Points:1" "Points:2" 8330 ./outputs/depth0_calvingVariables.txt 14  It appears that the velocity variables are given in x, y, z but I would like to calculate the magnitude. The next step is therefore to calculate the magnitude of velocity in 2D (just taking into account velocity in the x and y directions). [14]:  import numpy as np import sys # Get the header line and add 'velocity_magnitude', then save this as a file fp=open('./outputs/depth0_calvingVariables.txt') header = fp.readline().strip().split(' ') header.append("velocity_magnitude") np.savetxt('./outputs/calvingData_names.txt', header, delimiter=",", fmt="%s") # Calculate the velocity magnitude and append to dataset dataset = np.genfromtxt('./outputs/depth0_calvingVariables.txt', skip_header=1) temp=[] for row in dataset: velocity_mag = np.sqrt((float(row[3]))**2+(float(row[4])**2)) temp.append(velocity_mag) temp_array = np.array(temp).reshape(8329,1) # Save the new datatset result = np.append(dataset, temp_array, axis=1) np.savetxt('./outputs/calving_velocity_mag.csv', result, delimiter=",", fmt="%s")  In order to double check that the operation was successful, the number of columns will be counted and the first few lines of data will be printed. [15]:  %%bash tail -n 1 ./outputs/depth0_calvingVariables.txt | grep -o " " | wc -l tail -n 1 ./outputs/calving_velocity_mag.csv | grep -o "," | wc -l head ./outputs/calving_velocity_mag.csv  14 15 129.01,0.0,0.0,0.0,0.0,0.0,5.3266e-11,-5.3521e-11,-1.2158e-12,0.0,0.0,0.0,448350.0,8758700.0,7.1788,0.0 192.91,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,449270.0,8758200.0,6.9542,0.0 177.27,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,448800.0,8758500.0,10.856,0.0 148.28,0.0,0.0,-781.79,782.08,20.19,-2.3741,2.3852,0.054103,0.0,0.0,0.0,448300.0,8758600.0,-1.6698e-09,1105.823100907193 151.78,0.0,0.0,-1938.0,1134.2,-8.5292e-08,-5.922,3.506,1.2153e-10,0.0,0.0,0.0,448200.0,8758100.0,4.1405e-11,2245.4963014888267 162.41,0.0,0.0,-1073.8,993.68,-4.7606e-08,-3.2781,3.0476,6.7156e-11,0.0,0.0,0.0,448270.0,8758400.0,-2.3167e-12,1463.0264462408054 161.58,0.0,0.0,-1258.8,1170.4,-5.7489e-08,-3.8313,3.5924,7.3778e-11,0.0,0.0,0.0,448200.0,8758400.0,-2.2002e-11,1718.840772148485 152.99,0.0,0.0,-1741.2,1127.1,-7.164e-08,-5.3043,3.4774,1.2594e-10,0.0,0.0,0.0,448210.0,8758200.0,2.8983e-10,2074.1581063168737 155.49,0.0,0.0,-936.15,881.14,-5.2037e-08,-2.8785,2.6741,7.6945e-11,0.0,0.0,0.0,448290.0,8758500.0,4.1447e-11,1285.606674725983 154.17,0.0,0.0,-1629.5,1150.1,-7.2964e-08,-4.9626,3.5429,9.8018e-11,0.0,0.0,0.0,448200.0,8758200.0,-2.6549e-12,1994.492481810849  To get a coherent view of what the velocities look like across the glacier, I want to convert the point data into a raster. This can be done using gdal_grid but, in order to do this, a vrt header must be made for the data file. [16]:  %%bash cat ./inputs/calving_velocity_mag.vrt  <OGRVRTDataSource> <OGRVRTLayer name="calving_velocity_mag"> <SrcDataSource>CSV:./outputs/calving_velocity_mag.csv</SrcDataSource> <GeometryType>wkbPoint</GeometryType> <GeometryField encoding="PointFromColumns" x="field_13" y="field_14" z="field_16"/> </OGRVRTLayer> </OGRVRTDataSource>  The raster can now be created using gdal_grid. There are a number of different options available with regards to the interpolation algorithm. For this project, I have chosen to use the inverse distance to a power (‘invdist’) algorithm, which is a weighted average interpolator. Both the weighting power ‘p’ and the smoothing parameter ‘s’ can be user defined. As velocities across a glacier do not vary abruptly, I have opted for a high smoothing paramater of 500. The raster is then clipped by extent using a shapefile that covers the area of Kronebreen (via gdalwarp), before being visualised using openev. [17]:  %%bash gdal_grid -a invdist:power=1.0:smoothing=500.0 -outsize 400 400 -of GTiff -ot Float64 -l \ calving_velocity_mag ./inputs/calving_velocity_mag.vrt ./outputs/calving_velocity_mag.tiff gdalwarp -overwrite -cutline ./inputs/Domain.shp -crop_to_cutline -s_srs EPSG:32633 -t_srs EPSG:32633 \ ./outputs/calving_velocity_mag.tiff ./outputs/cropped_calving_velocity_mag.tiff -dstnodata 9999 /usr/bin/openev/bin/openev ./outputs/cropped_calving_velocity_mag.tiff  Grid data type is "Float64" Grid size = (400 400). Corner coordinates = (447926.775000 8764912.375000)-(466553.225000 8754987.625000). Grid cell size = (46.450000 24.750000). Source point count = 8329. Algorithm name: "invdist". Options are "power=1.000000:smoothing=500.000000:radius1=0.000000:radius2=0.000000:angle=0.000000:max_points=0:min_points=0:nodata=0.000000" 0...10...20...30...40...50...60...70...80...90...100 - done. Creating output file that is 522P x 279L. Processing ./outputs/calving_velocity_mag.tiff [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done. Default software rendering mode (use -h if accelerated video card installed). Loading tools from /usr/bin/openev/tools/Tool_DriverList.py Loading tools from /usr/bin/openev/tools/Tool_Export.py Loading tools from /usr/bin/openev/tools/Tool_ShapesGrid.py  These modelled velocities are now ready to be compared to the observed velocities, which are also provided as a .tiff. For this operation, I use gdal_calc to subtract one raster from the other and thus create a map of difference. I first use gdalwarp to make sure the two rasters have the same extent and resample the rasters so the pixel size is consistent. The extents are still slightly different, but this can be worked around by using the –extent flag in gdal_calc. It is notable that the two rasters have different units; observations are given in m/day whereas the modelled velocities are in m/year. The observed velocities are therefore multiplied by 365.25 within gdal_calc. [18]:  %%bash gdalwarp -overwrite -cutline ./inputs/Domain.shp -crop_to_cutline -s_srs EPSG:32633 -t_srs EPSG:32633 \ ./inputs/kroneVelo_Observed.tif ./outputs/kroneVelo_Observed_crop.tiff -dstnodata 9999 gdalwarp -tr 10 -10 ./outputs/kroneVelo_Observed_crop.tiff ./outputs/kroneVelo_Observed_crop_resample.tiff gdalwarp -tr 10 -10 ./outputs/cropped_calving_velocity_mag.tiff ./outputs/cropped_calving_velocity_mag_resample.tiff gdal_calc.py -A ./outputs/cropped_calving_velocity_mag_resample.tiff -B ./outputs/kroneVelo_Observed_crop_resample.tiff \ --outfile=./outputs/model_obs_comparison.tif --calc="A-(B*365.25)" --extent="intersect" /usr/bin/openev/bin/openev ./outputs/model_obs_comparison.tif  Creating output file that is 6309P x 3362L. Processing ./inputs/kroneVelo_Observed.tif [1/1] : 0Using internal nodata values (e.g. 0) for image ./inputs/kroneVelo_Observed.tif. ...10...20...30...40...50...60...70...80...90...100 - done. Creating output file that is 1863P x 993L. Processing ./outputs/kroneVelo_Observed_crop.tiff [1/1] : 0Using internal nodata values (e.g. 9999) for image ./outputs/kroneVelo_Observed_crop.tiff. Copying nodata values from source ./outputs/kroneVelo_Observed_crop.tiff to destination ./outputs/kroneVelo_Observed_crop_resample.tiff. ...10...20...30...40...50...60...70...80...90...100 - done. Creating output file that is 1859P x 993L. Processing ./outputs/cropped_calving_velocity_mag.tiff [1/1] : 0Using internal nodata values (e.g. 9999) for image ./outputs/cropped_calving_velocity_mag.tiff. Copying nodata values from source ./outputs/cropped_calving_velocity_mag.tiff to destination ./outputs/cropped_calving_velocity_mag_resample.tiff. ...10...20...30...40...50...60...70...80...90...100 - done. 0.. 0.. 1.. 2.. 3.. 4.. 5.. 5.. 6.. 7.. 8.. 9.. 10.. 10.. 11.. 12.. 13.. 14.. 15.. 15.. 16.. 17.. 18.. 19.. 20.. 20.. 21.. 22.. 23.. 24.. 25.. 25.. 26.. 27.. 28.. 29.. 30.. 30.. 31.. 32.. 33.. 34.. 35.. 35.. 36.. 37.. 38.. 39.. 40.. 40.. 41.. 42.. 43.. 44.. 45.. 45.. 46.. 47.. 48.. 49.. 50.. 50.. 51.. 52.. 53.. 54.. 55.. 55.. 56.. 57.. 58.. 59.. 60.. 60.. 61.. 62.. 63.. 64.. 65.. 65.. 66.. 67.. 68.. 69.. 70.. 70.. 71.. 72.. 73.. 74.. 75.. 75.. 76.. 77.. 78.. 79.. 80.. 80.. 81.. 82.. 83.. 84.. 85.. 85.. 86.. 87.. 88.. 89.. 90.. 90.. 91.. 92.. 93.. 94.. 95.. 95.. 96.. 97.. 98.. 99.. 100 - Done Default software rendering mode (use -h if accelerated video card installed). Loading tools from /usr/bin/openev/tools/Tool_DriverList.py Loading tools from /usr/bin/openev/tools/Tool_Export.py Loading tools from /usr/bin/openev/tools/Tool_ShapesGrid.py  [19]:  import rasterio from rasterio.plot import show from matplotlib import pyplot from matplotlib import colors, cm import numpy as np diff = rasterio.open("./outputs/model_obs_comparison.tif") modelled = rasterio.open("./outputs/cropped_calving_velocity_mag_resample.tiff") observed = rasterio.open("./outputs/kroneVelo_Observed_crop_resample.tiff") cmap = cm.bwr.copy() cmap.set_over('w',1.) cmap.set_under('w',1.) pyplot.imshow(diff.read(1), vmin=700, vmax=2500, cmap=cmap) pyplot.colorbar() pyplot.title('Modelled - Observed velocities (m/yr)') pyplot.show() pyplot.imshow(modelled.read(1), vmin=1500, vmax=2500, cmap=cmap) pyplot.colorbar() pyplot.title('Modelled velocities (m/yr)') pyplot.show() pyplot.imshow((observed.read(1))*365.25, vmin=0, vmax=1100, cmap=cmap) pyplot.colorbar() pyplot.title('Observed velocities (m/yr)') pyplot.show()  The final results from Part 1 show that modelled velocities are systematically higher than the observed velocities. This suggests that there may be an issue with a parameter such as basal friction or ice viscosity. It should be noted that some of the differences are explainable: 1. The very front of the glacier shows very low velocities. This is due to the retreat of the glacier and so this area should be ignored 2. Observed velocities contain error (est. up to 150 m/yr). This is particlarly true for slow flowing regions where little change occurs being subsequent satellite images. Therefore, we expect greater error over the upper parts of the glacier and place less emphasis on getting perfect correspondaance here. Despite this, the above analysis is of great help in terms of tuning and evaluating my ice sheet model set-up and has identified that I need to look into why the glacier is flowing faster than expected. ### 1.6.1.3. Part 2: Time series of Velocity and Surface Mass Balance I want to analyse a time series of frontal velocities from Kronebreen, Svalbard. I have configured the Elmer/Ice model to output the average velocity for each partition, with partition 0 corresponding to the frontal section of the glacier. As this is the most dynamic and interesting part of the glacier, it is this area that will be focused on for this project. I first want to explore the data files and then begin with a simple visualisation of how velocity varies throughout the simulation. The data files do not have headers, but a separate file ‘calving_averages.dat.names’ provides a list of the variables included. [20]:  %%bash # View the metadata - variable names, where the data comes from etc cat ./inputs/calving_averages/calving_averages.dat.names # Look at the first row of data head -n 1 ./inputs/calving_averages/calving_averages.dat.0 # Check the number of rows wc -l ./inputs/calving_averages/calving_averages.dat.0  Elmer version: 9.0 Elmer revision: ae1635a Elmer compilation date: 2021-01-20 Solver input file: krone_calving.sif File started at: 2021/05/19 17:31:54 Variables in columns of matrix: ./calving_averages.dat 1: mean: velocity 1 2: mean: velocity 2 3: mean: velocity 3 4: mean: smb 5: res: solver cpu time distance 6: res: solver real time distance 7: res: solver cpu time navier-stokes 8: res: solver real time navier-stokes 9: res: solver cpu time stresssolver 10: res: solver real time stresssolver 11: res: solver cpu time longitudinal mesh update 12: res: solver real time longitudinal mesh update 13: res: solver cpu time free surface top 14: res: solver real time free surface top 15: res: solver cpu time free surface bottom 16: res: solver real time free surface bottom 17: res: solver cpu time vertical mesh update 18: res: solver real time vertical mesh update 19: res: solver cpu time 3d calving 20: res: solver real time 3d calving 21: res: solver cpu time remesh 22: res: solver real time remesh 23: mean: velocity magnitude -1.06E+03,-1.30E+02,-2.49E+03,-3.12E+00,8.07E-03,1.58E-02,4.50E+01,4.63E+01,8.03E-01,8.27E-01,1.67E-01,1.81E-01,4.62E-01,5.49E-01,1.45E-01,1.76E-01,3.54E-01,3.89E-01,1.20E-01,4.07E-01,2.38E+00,3.08E+00,1064.157847 52 ./inputs/calving_averages/calving_averages.dat.0  [21]:  %%bash # Plot velocity magnitude R --vanilla -q # Prepare the data and create an array 'timesteps' #tseries = read.table("./inputs/longSpinUP_averages/longSpinUp_averages.dat.0", sep="," , header=FALSE) tseries = read.table("./inputs/calving_averages/calving_averages.dat.0", sep="," , header=FALSE) head(tseries) timesteps = t(1:52) #Create the plot png("./images/explore_velocities.png") plot(x=timesteps, y=tseries$V23, pch=4, cex=2, col="red", type='l', ylab="Velocity magnitude (m/yr)", xlab="Timestep")
dev.off()

#Spare code to plot x and y velocities...
#plot(x=timesteps, y=tseries$V3, pch=4, cex=2, col="blue", xlab=("Timestep"), ylab=("Y Velocity (blue) (m/yr)"), type='l', ylim=c(-550,1000)) #par(new=T) #plot(x=timesteps, y=tseries$V2, pch=4, cex=2, col="red", type='l', ylim=c(-3500,-600), axes=F, ylab=NA, xlab=NA)+axis(side=4)+mtext("X Velocity (red) (m/yr)", side=4,line=-1.0)

>
> # Prepare the data and create an array 'timesteps'
> #tseries = read.table("./inputs/longSpinUP_averages/longSpinUp_averages.dat.0", sep="," , header=FALSE)
> tseries = read.table("./inputs/calving_averages/calving_averages.dat.0", sep="," , header=FALSE)
V1     V2      V3      V4      V5      V6    V7    V8    V9   V10   V11
1 -1060 -130.0 -2490.0  -3.120 0.00807 0.01580 45.00 46.30 0.803 0.827 0.167
2 -3510  402.0   -53.7   5.930 0.00561 0.00937  5.84  6.08 0.172 0.186 0.108
3 -2160   -6.7 -2490.0   0.955 0.00351 0.00887  4.69  4.99 0.171 0.193 0.102
4 -3680  988.0   -66.2 -16.600 0.00569 0.00635  4.72  4.93 0.175 0.189 0.110
5 -3640  969.0   -70.2 -23.400 0.00503 0.00979  8.45  8.91 0.170 0.187 0.111
6 -3400  926.0   -59.0 -13.000 0.00435 0.01050  7.22  7.69 0.173 0.192 0.122
V12    V13    V14    V15    V16    V17    V18    V19   V20   V21  V22
1 0.181 0.4620 0.5490 0.1450 0.1760 0.3540 0.3890 0.1200 0.407 2.380 3.08
2 0.118 0.0349 0.0502 0.0201 0.0278 0.0675 0.0726 0.0868 0.216 0.713 1.04
3 0.114 0.0322 0.0507 0.0131 0.0186 0.0664 0.0720 0.0796 0.203 0.626 1.00
4 0.120 0.0383 0.0526 0.0272 0.0363 0.0660 0.0738 0.0764 0.206 0.720 1.76
5 0.122 0.0380 0.0544 0.0165 0.0279 0.0680 0.0738 0.0871 0.202 0.591 1.70
6 0.132 0.0423 0.0605 0.0157 0.0271 0.0702 0.0752 0.0876 0.221 0.781 1.09
V23
1 1064.158
2 3529.869
3 2163.036
4 3814.982
5 3767.336
6 3523.405
> timesteps = t(1:52)
>
> #Create the plot
> png("./images/explore_velocities.png")
> plot(x=timesteps, y=tseries$V23, pch=4, cex=2, col="red", type='l', ylab="Velocity magnitude (m/yr)", xlab="Timestep") > dev.off() null device 1 > > #Spare code to plot x and y velocities... > #plot(x=timesteps, y=tseries$V3, pch=4, cex=2, col="blue", xlab=("Timestep"), ylab=("Y Velocity (blue) (m/yr)"), type='l', ylim=c(-550,1000))
> #par(new=T)
> #plot(x=timesteps, y=tseries$V2, pch=4, cex=2, col="red", type='l', ylim=c(-3500,-600), axes=F, ylab=NA, xlab=NA)+axis(side=4)+mtext("X Velocity (red) (m/yr)", side=4,line=-1.0) >  [22]:  from IPython.display import Image Image("./images/explore_velocities.png")  [22]:  It looks like velocity begins high (c. 3500 m/yr) at the beginning of the simulation - probably due to the initial conditions imposed by the simulation. After around 5 days, the velocity drops abd begins oscillating around 2000 m/yr. The velocity varies considerably between timesteps, with lows of 500 m/yr and highs of around 3000 m/yr. This may be due to differences in surface mass balance (SMB) or in internal dynamics due to the impacts of, for example, calving events. Each timestep is one day, and the above plot covers a 52 day period corresponding the the first 52 days of the calendar year in 2016. The model is forced by a surface mass balance (SMB) dataset which is calculated as follows: SMB = PR - RU - SU - ER where RU = ME + RA - RF PR is total precipitation (including snowfall and rainfall), RU is meltwater runoff, SU is total sublimation and ER is erosion from drifting snow. Runoff is composed of liquid water from rain (RA) and melt (ME) that is not retained or refrozen in firn (RF). I want to look into the relationship between velocity and SMB by comparing the time series of velocities with the time series of SMB. The first step in this analysis is to plot the time series against eachother and compare them qualitatively. The SMB data is from Noël et al. (2020), with some pre-processing having occured before use during this course project. [23]:  %%bash # Plot X velocity alongside SMB R --vanilla -q # Prepare the velocity and SMB data and create an array 'timesteps' input_series = read.table("./inputs/calving_averages/calving_averages.dat.0", sep="," , header=FALSE) timesteps = t(1:52) #Create the plots png("./images/velocities_SMB.png") plot(x=timesteps, y=input_series$V23, pch=4, cex=2, col="red", type='l', ylim=c(0,4000), ylab="Velocity magnitude(red) (m/yr)", xlab="Timestep")
par(new=T)
plot(x=timesteps, y=input_series$V5, pch=4, cex=2, col="green", type='l', axes=F, ylab=NA, xlab=NA)+axis(side=4)+mtext("SMB (green) (mm w.e./yr)", side=4,line=-1.0) dev.off()  > > # Prepare the velocity and SMB data and create an array 'timesteps' > input_series = read.table("./inputs/calving_averages/calving_averages.dat.0", sep="," , header=FALSE) > timesteps = t(1:52) > > #Create the plots > png("./images/velocities_SMB.png") > plot(x=timesteps, y=input_series$V23, pch=4, cex=2, col="red", type='l', ylim=c(0,4000), ylab="Velocity magnitude(red) (m/yr)", xlab="Timestep")
> par(new=T)
> plot(x=timesteps, y=input_series\$V5, pch=4, cex=2, col="green", type='l', axes=F, ylab=NA, xlab=NA)+axis(side=4)+mtext("SMB (green) (mm w.e./yr)", side=4,line=-1.0)
numeric(0)
> dev.off()
null device
1
>

[24]:

from IPython.display import Image
Image("./images/velocities_SMB.png")

[24]:

I want to see if there is a correlation between these two time series, and will start with a simple linear regression:

[25]:

%%bash
# Regress SMB against X velocity

R --vanilla -q

regression = lm(V23 ~ V5, data = input_series)

summary(regression)

>
> input_series = read.table("./inputs/calving_averages/calving_averages.dat.0", sep="," , header=FALSE)
>
> regression = lm(V23 ~ V5, data = input_series)
>
> summary(regression)

Call:
lm(formula = V23 ~ V5, data = input_series)

Residuals:
Min       1Q   Median       3Q      Max
-1457.71  -497.13   -53.78   452.40  1772.44

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   2372.8      352.5   6.732 1.58e-08 ***
V5          -58044.9    67071.9  -0.865    0.391
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 764.8 on 50 degrees of freedom
Multiple R-squared:  0.01476,   Adjusted R-squared:  -0.004947
F-statistic: 0.7489 on 1 and 50 DF,  p-value: 0.3909

>

The p value is poor (0.39) and shows the result is not statistically significant. The R^2 value is low at 0.15 suggesting that there is little correlation between surface mass balance and velocity magnitude. However, this is a very simply done analysis and is conducted only as an initial data exploration.

I will now try dynamic time warping as an alternative method for comparing the two time series:

[26]:

import numpy as np
from dtw import *
import matplotlib.pyplot as plt

# Get the data

query = np.loadtxt(fname="./inputs/calving_averages/calving_averages.dat.0", delimiter=',', usecols=22)
# Smb in mm w.e. per YEAR converted to ice equiv. (as used in Elmer/Ice)
reference = (np.loadtxt(fname="./inputs/calving_averages/calving_averages.dat.0", delimiter=',', usecols=4))*365.25/(1000.0/910.0)

# Perform the analysis

alignment = dtw(query, reference, keep_internals=True)

## Display the warping curve, i.e. the alignment curve
alignment.plot(type="threeway")
plt.show()

dtw(query, reference, keep_internals=True, window_type="sakoechiba", window_args={"window_size":10},
step_pattern=rabinerJuangStepPattern(6, "c")).plot(type="twoway",offset=-1)

alignment.normalizedDistance

[26]:

1374.9222062850727

The patterns of the two time series looked similar when plotted side by side but the distance between them is relatively large, which could be partially due to the very different units.

A very close correspondance is not expected as ice dynamics do not react to external changes immediately but are partly controlled by preceding conditions.

However, noise may still be impacting the result, and the warp path jumps (rather than being broadly diagonal) suggesting a less-than-optimal fit. In order to try and remove some noise, I will apply a wavelet transform to extract the key features of the two time series and then re-run the dynamic time warp.

[27]:

import pywt
# Try out wavelet feature extraction ...

#Velocity

sLst = np.arange(1, 31)
wav_transform_velocity, freqs = pywt.cwt(query, sLst, 'mexh')

plt.plot(wav_transform_velocity[1])
plt.plot(query)
plt.show()

#SMB

wav_transform_smb, freqs = pywt.cwt(reference, sLst, 'mexh')

plt.plot(wav_transform_smb[1])
plt.plot(reference)
plt.show()

[28]:

import numpy as np
from dtw import *
import matplotlib.pyplot as plt

# Run dynamic time warp again

alignment = dtw(wav_transform_velocity[1], wav_transform_smb[1], keep_internals=True)

## Display the warping curve, i.e. the alignment curve
alignment.plot(type="threeway")
plt.show()

dtw(wav_transform_velocity[1], wav_transform_smb[1], keep_internals=True, window_type="sakoechiba", window_args={"window_size":10},
step_pattern=rabinerJuangStepPattern(6, "c")).plot(type="twoway",offset=-1)

alignment.normalizedDistance

[28]:

468.9828699904016

With the wavelet pre-processing, there is a much reduced distance between the two time series (1375 vs 469), suggesting that there may some connection between the two. However, the warp path still does’nt take a consistent, diagonal path. This means that the analysis is not overly successful and other techniques may be required to investigate the relationship.

The lack of good correspondance points to the fact that other variables besides SMB are also influencing the velocity of the glacier. These could be indirect effects of SMB (for example through thinning and thickening of the glacier in different areas altering the gravitational driving stress), or unrelated variables such as bed geometry.

It is also notable that SMB variations are not expected to immediately cause large velocity variations, as mentioned before. However, understanding how quickly the glacier can respond to external perturbations is of interest with regards to reducing uncertainty regarding how they will behave in the future and so in improving sea level rise estimates.

The results of the qualitative and dynamic time warp analyses show that there may be a connection between the two, although nothing significant has been identified yet.

I next want to split up SMB into its underlying components, so I can use these time series in my analysis. This will be useful for understanding whether frontal velocities are more sensitive to certain parameters and is also useful in investigating whether a single/few parameter(s) could be used to predict dynamical changes.

The variables I will focus on are: 1) Precipitation (mm w.e.) 2) Air Temperature (K) 3) Snowfall (mm w.e.) 4) Snowmelt (mm w.e.

The variables are given in separate .tiff files which each contain 365 bands with the mean value for each day in the period 2000-2019. Data will be extracted from each file for relevant location(s) using gdallocationinfo.

Due to the poor results obtained by the dynamic time warp analysis, I want to instead move forward by using the Random Forest algorithm to try and fit a regression model on the data. In this way, I can understand the relative importance of the various factors for frontal velocity at Kronebreen.

[29]:

from osgeo import gdal
import os
import numpy as np
import sys

total=0.0
with open("./outputs/smb_components.csv", 'w') as myfile:
for band in range(1, 53):
for var in 'precip', 't2m', 'snowfall', 'snowmelt':
for x in range(448209, 449110, 110):
for y in range(8755649, 8758086, 250):
temp=os.popen('gdallocationinfo -valonly -geoloc -b %s ./inputs/daily_%s_means_2000to2019_CROPPED.tif %s %s' % (band, var, x, y)).read()
total=float(temp)+total
average=total/100
myfile.write(str(average)+', ')
total=0.0
myfile.write('\n')
myfile.close()

[30]:

import pandas as pd
import numpy as np
from sklearn.ensemble import RandomForestRegressor as RFReg
from sklearn.model_selection import train_test_split,GridSearchCV
from sklearn.pipeline import Pipeline
from scipy.stats.stats import pearsonr
import matplotlib.pyplot as plt

[31]:

vel_dataset = pd.read_csv("./inputs/calving_averages/calving_averages.dat.0", header=None)

#SMB components
x = smb_dataset.iloc[:, 0:4].values
# Velocity magnitude
y = vel_dataset.iloc[:, 22:23].values

variables = ['precip', 't2m', 'snowfall', 'snowmelt']

[32]:

x.shape

[32]:

(52, 4)

[33]:

y.shape

[33]:

(52, 1)

The data is split into training and testing datasets (various splits). After this, a pipeline is defined. This allows many different parameter combinations and scorings to be tested out.

After initial tests (not shown here), I added lower values for max_features and max_depth as well as a higher number for n_estimators into my pipeline.

In addition, I played around with different proportions of data used for training (20%, 50%, 70%) and different ways of scoring the regressions (R^2, explained variance, and negative root mean square error). The results from these tests are shown below:

Test size | Scoring | Train score (pearson) | Test score (pearson)|
— | — | — | |
0.2 | R^2 | 0.820 | -0.043 |
0.5 | R^2 | 0.916 | -0.094 |
0.7 | R^2 | 0.919 | 0.052 |
0.2 | Explained variance | 0.922 | -0.087 |
0.5 | Explained variance | 0.918 | -0.017 |
0.7 | Explained variance | 0.931 | 0.095 |
0.2 | Neg root mean square error | 0.816 | -0.068 |
0.5 | Neg root mean square error | 0.937 | -0.079 |
0.7 | Neg root mean square error | 0.951 | 0.091 |

The parameter combination which gives the best fit will then be used for analysis. This corresponds to: 1) A test size of 0.7 2) Scoring based on negative root mean square error 3) Max depth: 25 4) Max features: 0.33 5) Max samples: 0.7 6) Number of estimators: 500

[34]:

# Split into testing and training data sets
X_train, X_test, Y_train, Y_test = train_test_split(x, y, test_size=0.7, random_state=24)
y_train = np.ravel(Y_train)
y_test = np.ravel(Y_test)

[35]:

# Determine a pipeline to test out different parameters and find the set of
# parameters with the best fit. Various parameters are used to score the different fits.

pipeline = Pipeline([('rf',RFReg())])

parameters = {
'rf__max_features':("log2","sqrt",0.33, 0.2),
'rf__max_samples':(0.3, 0.5,0.6,0.7),
'rf__n_estimators':(500,1000,2000,3000),
'rf__max_depth':(25,50,100,200)}

grid_search = GridSearchCV(pipeline,parameters,n_jobs=-1,cv=3,scoring='neg_root_mean_squared_error',verbose=1)
grid_search.fit(X_train,y_train)

Fitting 3 folds for each of 256 candidates, totalling 768 fits

[35]:

GridSearchCV(cv=3, estimator=Pipeline(steps=[('rf', RandomForestRegressor())]),
n_jobs=-1,
param_grid={'rf__max_depth': (25, 50, 100, 200),
'rf__max_features': ('log2', 'sqrt', 0.33, 0.2),
'rf__max_samples': (0.3, 0.5, 0.6, 0.7),
'rf__n_estimators': (500, 1000, 2000, 3000)},
scoring='neg_root_mean_squared_error', verbose=1)

[36]:

grid_search.best_score_

[36]:

-859.084465363649

[37]:

# The best parameters are printed, and these are then used for the regression

print ('Best Training score: %0.3f' % grid_search.best_score_)
print ('Optimal parameters:')
best_par = grid_search.best_estimator_.get_params()
for par_name in sorted(parameters.keys()):
print ('\t%s: %r' % (par_name, best_par[par_name]))

Best Training score: -859.084
Optimal parameters:
rf__max_depth: 25
rf__max_features: 0.33
rf__max_samples: 0.7
rf__n_estimators: 500

[39]:

rfReg = RFReg(n_estimators=500,max_features=0.33,max_depth=25,max_samples=0.7,n_jobs=-1,random_state=24)
rfReg.fit(X_train, y_train);
dic_pred = {}
dic_pred['train'] = rfReg.predict(X_train)
dic_pred['test'] = rfReg.predict(X_test)
[pearsonr(dic_pred['train'],y_train)[0],pearsonr(dic_pred['test'],y_test)[0]]

[39]:

[0.950879197769439, 0.09090151937598724]

Using Pearson’s R, it can be seen that there is agood fit for the testing data (0.951) but a very poor fit for the testing data (0.091).

This is likely due to overfitting and the following parameters could be altered to try and fix this:

1. max_features - Reducing this number will reduce the number of features each tree is randomly assigned

2. max_depth - Reducing this parameter lowers the complexity of the learned models, and thus minimises the risk of over fitting

3. n_estimators - Overfitting is less likely with more estimators

However,despite adding more possible parameters into my pipeline, no solution could be found to improve the results and prevent overfitting. A key limitation for this model is the lack of input data. Although this is a dissapointing realisation, it does give me hope that when I move on to running longer simulations (with e.g. over 1000 timesteps), I may see some progress.

Despite this, I still want to visualise the results and look into the relative importance of all the various predictors for the best fit regression that I have:

[40]:

plt.scatter(y_train,dic_pred['train'])
plt.xlabel('orig')
plt.ylabel('pred')
ident = [0, 4000]
plt.plot(ident,ident,'r--')

[40]:

[<matplotlib.lines.Line2D at 0x7f0a5200e940>]

[41]:

plt.scatter(y_test,dic_pred['test'])
plt.xlabel('orig')
plt.ylabel('pred')
ident = [0, 4000]
plt.plot(ident,ident,'r--')

[41]:

[<matplotlib.lines.Line2D at 0x7f0a5202ec10>]

[42]:

# Get numerical feature importances
importances = list(rfReg.feature_importances_)
# List of tuples with variable and importance
feature_importances = [(feature, round(importance, 2)) for feature, importance in zip(variables, importances)]
# Sort the feature importances by most important first
feature_importances = sorted(feature_importances, key = lambda x: x[1], reverse = True)
# Print out the feature and importances
[print('Variable: {:20} Importance: {}'.format(*pair)) for pair in feature_importances];

Variable: precip               Importance: 0.28
Variable: t2m                  Importance: 0.25
Variable: snowfall             Importance: 0.24
Variable: snowmelt             Importance: 0.23

[45]:

# Plot the importance of different predictors
impt = [rfReg.feature_importances_, np.std([tree.feature_importances_ for tree in rfReg.estimators_],axis=1)]
ind = np.argsort(impt[0])
plt.rcParams["figure.figsize"] = (10,6)
plt.barh(range(len(variables)),impt[0][ind],color="r", xerr=impt[1][ind], align="center")
plt.yticks([3, 2, 1, 0], ['precipitation', 'temperature', 'Snowfall', 'Snowmelt'])

[45]:

([<matplotlib.axis.YTick at 0x7f0a51e69220>,
<matplotlib.axis.YTick at 0x7f0a51e78b80>,
<matplotlib.axis.YTick at 0x7f0a51e67a30>,
<matplotlib.axis.YTick at 0x7f0a51de0a00>],
[Text(0, 3, 'precipitation'),
Text(0, 2, 'temperature'),
Text(0, 1, 'Snowfall'),
Text(0, 0, 'Snowmelt')])

The final results show that the predictors used in this study are not sufficient to predict frontal velocity at Kronebreen (at least for this short time period). However, all the predictors did have some importance and so are not irrelevant parameters. The relative importance of the different parameters was similar, but snowmelt came out as being the most important predictor.

Full reference for SMB data:

Noël, Brice P Y; Jakobs, Constantijn L; van Pelt, Ward; Lhermitte, Stef; Wouters, Bert; Kohler, Jack; Hagen, Jon Ove; Luks, Bartłomiej; Reijmer, Carleen; van de Berg, Willem Jan; van den Broeke, Michiel R (2020): Annual surface mass balance (SMB) and components of Svalbard glaciers statistically downscaled to 500 m spatial resolution (1958-2018). PANGAEA, https://doi.org/10.1594/PANGAEA.920984


text:
# Functional Analysis II Spring 2018

Lecturer
Prof. Dr. Alessandro Carlotto
Course Assistant
Mario Schulz
Teaching Assistants
Alessandro Pigati
Xavier Fernandez-Real Girona
Lectures
Monday 10–12 / HG G 5
Thursday 13–15 / HG G 5
Exercise classes
Monday 09–10
Office hours
Tuesday 14–15
Wednesday 13–15
First Lecture
19.02.2018
First Exercise class
19.02.2018
Course Catalogue
401-3462-00L Functional Analysis II
Prerequisite
Functional Analysis I

### Primary references

Michael Struwe. Funktionalanalysis I und II. Lecture notes, ETH Zürich, 2013/14.

Haim Brezis. Functional analysis, Sobolev spaces and partial differential equations. Universitext. Springer, New York, 2011.

### Extra references

David Gilbarg, Neil Trudinger. Elliptic partial differential equations of second order. Classics in Mathematics. Springer-Verlag, Berlin, 2001.

Qing Han, Fanghua Lin. Elliptic partial differential equations. Second edition. Courant Lecture Notes in Mathematics, 1. Courant Institute of Mathematical Sciences, New York; American Mathematical Society, Providence, RI, 2011.

Michael Taylor. Partial differential equations I. Basic theory. Second edition. Applied Mathematical Sciences, 115. Springer, New York, 2011.

Lars Hörmander. The analysis of linear partial differential operators. I. Distribution theory and Fourier analysis. Classics in Mathematics. Springer-Verlag, Berlin, 2003.

DateContentReference 1 (Struwe)Reference 2 (Brezis)Extras
19.02.2018 Introduction of the course. A model problem: the elastic membrane with fixed boundary.
The Euler-Lagrange equation associated to a functional. A general roadmap to elliptic regularity.
–  –  –
22.02.2018 Solvability of the Poisson equation: the Poincaré inequality and two existence results for weak solutions,
either via Riesz theorem or by Fredholm alternative.
§ 7.1 § 8.1,  § 8.3  –
26.02.2018 Distributional and weak derivatives, examples and basic facts.
Definition of Sobolev spaces. Completeness, separability and reflexivity.
§ 7.2 § 9.1  –
01.03.2018 Absolute continuity of functions in $$W^{1,p}(I)$$, weak vs. pointwise derivative.
(Dual) characterizations of $$W^{1,p}(I)$$ for $$p>1$$. The Sobolev extension operator.
§ 7.3
(till Satz 7.3.4)
§ 8.2
(till Thm 8.6)
–
05.03.2018 Density of test functions in $$W^{1,p}(I)$$ for all $$1\leq p < \infty$$.
The one-dimensional Sobolev embedding theorem. $$W^{1,p}$$ as a Banach algebra.
§ 7.3 § 8.2  –
08.03.2018 Discussion of some examples of second-order ODEs with either Dirichlet or Neumann boundary conditions.
Perturbation methods and convergence of approximation schemes.
§ 7.4 § 8.4  –
12.03.2018 Sobolev functions of $$N$$ variables, examples of singular functions in $$W^{1,p}$$. A criterion for the coincidence of pointwise and weak derivative (null capacity). The Meyers-Serrin approximation theorem, and comments on boundary behaviour. § 8.1,  § 8.2 § 9.1  –
15.03.2018 Equivalent characterizations of $$W^{1,p}$$. Lipschitz vs. $$W^{1,\infty}$$.
Calculus rules for Sobolev functions: sums, products and compositions (chain rule).
§ 8.3 § 9.1  –
19.03.2018 Extension of functions in $$W^{1,p}(\Omega)$$ for $$\Omega$$ a relatively compact domain of class $$C^1$$.
Two approximation theorems for $$W^{1,p}$$ functions, on bounded or unbounded domains.
§ 8.4
(till Satz 8.4.2)
§ 9.2  –
22.03.2018 Imposing boundary values for elliptic problems: trace operators and their properties.
Two characterizations of $$H^1_0(\Omega)$$ and the canonical splitting of $$H^1(\Omega)$$.
§ 8.4 complements to § 9  –
26.03.2018 The Sobolev-Gagliardo-Nirenberg inequality, link with the isoperimetric inequality in $$\mathbb{R}^n$$.
The Sobolev embedding theorem for $$p < n$$, comments on the borderline case $$p = n$$.
§ 8.5 § 9.3  –
29.03.2018 A panorama on the Sobolev embedding theorems, exponential integrability of elements in $$W^{1,n}(\Omega)$$. Hölder continuity, completeness of the associated functional spaces, compact embedding of $$C^{0,\alpha}(\Omega)\subset C^{0,\beta}(\Omega)$$ for $$\alpha>\beta$$. § 8.6.1  –  –
09.04.2018 Morrey-Campanato spaces, integral characterization of Hölder continuity.
The Sobolev embedding theorem for $$p>n$$ and associated compactness results.
§ 8.6.2 § 9.3  –
12.04.2018 Pointwise differentiability of functions in $$W^{1,p}(\mathbb{R}^n)$$ for any $$p>n$$.
Higher-order Sobolev embedding theorems and their applications to regularity of weak solutions.
§ 8.6.3-4,  § 9.1 § 9.3,  § 9.5  –
16.04.2018 Interior regularity for solutions of the Poisson equation: formal $$H^1$$ and $$H^2$$ estimates, and their rigorous counterpart via Nirenberg's method (difference quotients). § 9.2
except § 9.2.1
§ 9.6  –
19.04.2018 Higher Sobolev estimates for weak solutions of the Poisson equation.
Ellipticity for operators in divergence form, and corresponding $$H^{k+2}$$ interior estimates.
§ 9.2,  § 9.5 § 9.6  –
23.04.2018 Sobolev estimates and boundary regularity for weak solutions of the Poisson equation with flat boundary.
The case of curved boundary: flattening via diffeomorphisms and the modified equation.
§ 9.3
(till Lem 9.3.4)
§ 9.6  –
26.04.2018 Sobolev estimates and boundary regularity for weak solutions of elliptic equations.
Global $$H^{k+2}$$ estimates for weak solutions of elliptic equations.
§ 9.3 § 9.6  –
30.04.2018 The closure of $$-\Delta\colon C^2_0\subset L^2\to L^2$$.
The Dirichlet spectrum of the Laplace operator and the min-max characterization of its eigenvalues.
§ 9.4 § 9.8  –
03.05.2018 Introduction to Schauder theory, context and motivations.
Campanato estimates for solutions of homogeneous problems (with constant coefficients).
§ 10.1,  § 10.2  –  –
07.05.2018 Campanato estimates for solutions of homogeneous problems (with constant coefficients): the boundary case.
Estimates for solutions of inhomogeneous problems (with constant coefficients).
§ 10.2  –  –
14.05.2018 Morrey spaces and their equivalence to Campanato spaces for $$\nu < n$$.
Local $$C^{1,\alpha}$$ Schauder estimates (interior and boundary cases).
§ 10.3,  § 10.4
(till Lem 10.4.3)
–  –
17.05.2018 Local $$C^{2,\alpha}$$ Schauder estimates (interior and boundary cases).
Global $$C^{2,\alpha}$$ Schauder estimates.
§ 10.4  –  –
24.05.2018 Solvability in $$C^{2,\alpha}$$ for elliptic boundary value problems.
The case of the Laplacian and the method of continuity.
§ 10.5  –  –
28.05.2018 The weak maximum principle for elliptic operators, and applications to the inhomogeneous case.
Connections with the maximum modulus principle and the Schwarz Lemma for holomorphic functions.
–  – weak maximum principle
31.05.2018 The strong maximum principle via the Hopf boundary point lemma. Elliptic barriers.
The method of sub- and super-solutions, examples of geometric relevance.
–  – strong maximum principle
Dates
6–9 August 2018
Place
HG G 26.1
 Rules of the exam: Rules.pdf Advice for the exam: Advice.pdf Sample Questions: Samplequestions.pdf

Please register in myStudies and use echo.ethz.ch to enroll for an exercise class.

AssistantClassroom
Xavier Fernandez-Real Girona HG G 26.3
Alessandro Pigati HG E 33.3
Mario Schulz HG F 26.5

Every Thursday, a new problem set is uploaded here. You have one week time to solve the problems. The following Monday during exercise class the assistants show the right approaches for the more involved problems. If you still have difficulties understanding or solving certain tasks, please prepare your questions and join the office hours on Tuesday or Wednesday.

#### Homework collection and delivery:

You may hand in your homework for grading to the assistant of the exercise class you enrolled for in myStudies. Either leave it in the respective letterbox in HG F28 or – if you wrote your homework in LaTeX – you may send the pdf-file to your assistant by email. The deadline is on Thursdays 4 pm. The graded homework sheets are delivered in class on Thursdays by the instructor.

Every problem is marked by one of the following symbols.

Computation
Get your hands dirty and calculate.
Bookkeeping
Apply what you learn in basic situations.
Comprehension
Construct examples and give full proofs.
Hard problem
Challenging problems are denoted by one up to three diamonds. It is recommended that you start working on these problems only after you have reviewed the weekly material and carefully solved all other exercises in the assignment.
19.02.2018 review
22.02.2018 01.03.2018 Problem Set 1 Hints 1 Solution 1 updated 28.02.18
01.03.2018 08.03.2018 Problem Set 2 Hints 2 Solution 2
08.03.2018 15.03.2018 Problem Set 3 Hints 3 Solution 3
15.03.2018 22.03.2018 Problem Set 4 Solution 4 Online quiz: echo.ethz.ch
22.03.2018 29.03.2018 Problem Set 5 Hints 5 Solution 5
29.03.2018 12.04.2018 Problem Set 6 Hints 6 Solution 6
12.04.2018 19.04.2018 Problem Set 7 Hints 7 Solution 7
19.04.2018 26.04.2018 Problem Set 8 Solution 8 Online quiz: echo.ethz.ch update: 8.12. convergence of a subsequence
26.04.2018 03.05.2018 Problem Set 9 Hints 9 Solution 9 updated 30.04.18
03.05.2018 11.05.2018 Problem Set 10 Solution 10 Online quiz: echo.ethz.ch
11.05.2018 24.05.2018 Problem Set 11 Hints 11 Solution 11
25.05.2018 31.05.2018 Problem Set 12 Solution 12 Online quiz: echo.ethz.ch

You are free to come and ask questions. The location changes depending on which assistant is on duty. The schedule is as follows (up to possible short-term changes, please check for updates):

DayDateTimeLocationAssistant
Tuesday20.02.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday21.09.201813:00–15:00 HG G 28 Xavier Fernandez-Real Girona
Tuesday27.02.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday28.02.201813:00–15:00 HG F 28.3 Mario Schulz
Tuesday06.03.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday07.03.201813:00–15:00 HG G 28 Xavier Fernandez-Real Girona
Tuesday13.03.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday14.03.201813:00–15:00 HG F 28.3 Mario Schulz
Tuesday20.03.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday21.03.201813:00–15:00 HG G 28 Xavier Fernandez-Real Girona
Tuesday27.03.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday28.03.201813:00–15:00 HG F 28.3 Mario Schulz
Tuesday10.04.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday11.04.201813:00–15:00 HG G 28 Xavier Fernandez-Real Girona
Tuesday17.04.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday18.04.201813:00–15:00 HG F 28.3 Mario Schulz
Tuesday24.04.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday25.04.201813:00–15:00 HG G 28 Xavier Fernandez-Real Girona
Wednesday02.05.201813:00–15:00 HG F 28.3 Mario Schulz
Tuesday08.05.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday09.05.201813:00–15:00 HG G 28 Xavier Fernandez-Real Girona
Tuesday15.05.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday16.05.201813:00–15:00 HG F 28.3 Mario Schulz
Tuesday22.05.201814:00–15:00 HG FO 27.9 Alessandro Pigati
Wednesday23.05.201813:00–15:00 HG G 28 Xavier Fernandez-Real Girona
Tuesday29.05.201814:00–15:00 HG F 28.3 (update) Mario Schulz
Wednesday30.05.201813:00–15:00 HG F 28.3 Mario Schulz


text:
# Model Evaluation: a crucial step in solving a machine learning problem

#### Machine Learning (ML) model evaluation

Machine Learning has developed several efficient techniques for a wide range of problems. A single technique or model like Googlenet is used across various problems such as image classification, object detection and others. Some models like MobileNet are designed to work efficiently in computational limited resources like mobile devices. It is a challenge to figure out the technique/ model that will be the best fit for your problem. To solve this, we will show how to evaluate models to determine if its performance meets your expectations/ requirements.

There are 3 major steps in solving a machine learning problem:

• Apply machine learning algorithms to solve the problem

When you reach step 3, a major issue arises: Which algorithm is best suited to solve our problem?

There are a large number of machine learning algorithms out there but not all of them apply to a given problem. We need to choose among those algorithms the one that best suits our problem and gives us the desired results. This is where the role of Model Evaluation comes in. It defines metrics to evaluate our models and then based on that evaluation, we choose one or more than one model to use. Let’s see how we do it.

We evaluate a model based on:

• Test Harness
• Performance Measure
• Cross-validation
• Testing Algorithms

### Test Harness

The test harness is the data on which you will train and test your model against a performance measure. It is crucial to define which part of the data will be used for training the model and which part for testing it. This may be as simple as selecting a random split of data (66% for training, 34% for testing) or may involve more complicated sampling methods.

While training the model on the training dataset, it is not exposed to test dataset. Its predictions on the test dataset are indicative of the performance of the model in general.

### Performance Measure

The performance measure is the way you want to evaluate a solution to the problem. It is the measurement you will make of the predictions made by a trained model on the test dataset.

There are many metrics used for performance measure according to the type problem that we’re solving. Let’s see some of them:

• Classification Accuracy: It is the ratio of the number of correct predictions to the total number of input samples.
$$Accuracy = \frac{Number\ of\ correct\ predictions}{Total\ number\ of\ predictions\ made}$$
This metric is used only in those cases where we have an equal or nearly equal number of data points belonging to all the classes.

For example, if we have a dataset where 98 percent of the dataset belongs to class A and 2 percent to class be, then a model which simply predicts a given data point to be belonging to class A all the time will have 98 percent accuracy. But in reality, such a model is performing very poorly.

• Logarithmic Loss: Suppose, there are N samples belonging to M classes, then the Log Loss is calculated as below :
$$Logarithmic Loss = -\frac{1}{N}\sum_{i=0}^N\sum_{j=0}^M\ {y_i}_j*log({p_i}_j)$$

where,
${y_i}_j$ indicates whether sample i belongs to class j or not
${p_i}_j$ indicates the probability of sample i belonging to class j
This function generates values between [0, $\infty$). Log loss nearer to zero indicates higher accuracy and log loss away from zero indicates lower accuracy. It is mostly used in multiclass classification and gives fairly good results.

• Confusion Matrix: Let's assume we have a binary classification problem. We have some samples belonging to two classes: YES or NO. Also, we have our own classifier which predicts a class for a given input sample. On testing our model on 160 samples, we get the following result.

There are 4 important terms :

• True Positives: The cases in which we predicted YES and the actual output was also YES.
• True Negatives: The cases in which we predicted NO and the actual output was NO.
• False Positives: The cases in which we predicted YES and the actual output was NO.
• False Negatives: The cases in which we predicted NO and the actual output was YES.

Accuracy for the matrix can be calculated by taking the average of the values lying across the “main diagonal” i.e
$$Accuracy = \frac{True\ Positives\ +\ True Negatives}{Total\ number\ of\ samples}$$
i.e., $$Accuracy = \frac{40+70}{165}$$

• Mean Squared Error: Mean Squared Error(MSE) is quite similar to Mean Absolute Error, the only difference being that MSE takes the average of the square of the difference between the original values and the predicted values. The advantage of MSE is that it is easier to compute the gradient, whereas Mean Absolute Error requires complicated linear programming tools to compute the gradient. As we take the square of the error, the effect of larger errors become more pronounced then smaller error, hence the model can now focus more on the larger errors.

$$Mean Squared Error = \frac{1}{N} \sum_{j=1}^N (y_j\ -\ \hat y_j)^2$$

### Cross-validation

Here we use the entire dataset to train the model and test the model as well. Here’s how.

Step 1: we divide our dataset into equally sized groups of data points called folds.
Step 2: Then we train our data on all the folds except 1.
Step 3: Next we test our data on that fold that was left out.
Step 4: Repeat step 2 and 3 such that all the folds get to be the test data one by one.
Step 5: The performance measures are averaged across all folds to estimate the capability of the algorithm on the problem.

For example, if we have 3 folds.

Step 1: Train on fold 1+2, test on fold 3
Step 2: Train on fold 2+3, test on fold 1
Step 3: Train on fold 3+1, test on fold 2
Step 4: Average the results of all the above steps to get the final result.

Usually, the number of folds we use are 3, 5, 7 and 10 folds.

### Testing Algorithms

Now, its time to test a variety of machine learning algorithms.

Select 5 to 10 standard algorithms that are appropriate for your problem and run them through your test harness. By standard algorithms, one means popular methods no special configurations. Appropriate for your problem means that the algorithms can handle regression if you have a regression problem.

If you have quite a few numbers of algorithms to test, then you may want to reduce the size of the dataset.

Following the above steps, you will be able to determine the best algorithm for your system configuration and problem.


text:
MCQ Questions for Class 7 Maths Chapter 12 Algebraic Expressions with Answers

We have compiled the NCERT MCQ Questions for Class 7 Maths Chapter 12 Algebraic Expressions with Answers Pdf free download covering the entire syllabus. Practice MCQ Questions for Class 7 Maths with Answers on a daily basis and score well in exams. Refer to the Algebraic Expressions Class 7 MCQs Questions with Answers here along with a detailed explanation.

Algebraic Expressions Class 7 MCQs Questions with Answers

Choose the correct option.

Question 1.
The coefficient of x² in – 5x²y is
(a) -5
(b) y
(c) 5y
(d) -5y

Question 2.
The expression 4a²b – 5ab² + 1 is a
(a) monomial
(b) binomial
(c) trinomial
(d) None of these

Question 3.
What is the coefficient of x in the expression 5x³ + 7x² – 5?
(a) 5
(b) 7
(c) -5
(d) 0

Question 4.
Which of the following pairs of terms is of unlike terms
(a) -p²q², 12p²q²
(b) – 4x²; – Axy²
(c) 41, 100
(d) qp², 13p²q

Answer: (b) – 4x²; – Axy²

Question 5.
How many terms are there in the expression 5x + 3y + 8?
(a) 1
(b) 2
(c) 3
(d) 4

Question 6.
Subtracting 5x + 3y from 5x – y + 8 given
(a) 4x – 8
(b) 4y – 8
(c) 8 – 4y
(d) 4x + 8

Question 7.
The value of 3x² + 2x – 7 at x = -2 is
(a) 1
(b) -1
(c) 9
(d) 23

Question 8.
Subtracting -3x² – 1 from 2x + 1, we get
(a) – 3x² + 2x – 2
(b) – 3x² + 2
(c) 3x² + 2x + 2
(d) 3x² – 2

Answer: (c) 3x² + 2x + 2

Question 9.
The value of the expression
100 – 10x³ + 5x² at x = 0 is
(a) 10
(b) 100
(c) -10
(d) – 100

Question 10.
Simplify: – z² + 22z² – 5z – 11z² + 5z
(a) z²
(b) -z²
(c) 5z
(d) -5z

Question 11.
The coefficient of x² in -2x²y is:
(a) -2
(b) y
(c) 2y
(d) -2y

Question 12.
The expression 4p²q – 4pq² + 1 is a
(a) monomial
(b) binomial
(c) trinomial
(d) none of these

Fill in the blanks

Question 1.
The degree of the polynomial 5x³ – 3x² + 6x4 + 7 ……………….

Question 2.
An algebraic expression having only two terms is called a ……………..

Question 3.
The value of x + 6 when x = -5 is ………………

Question 4.
Constants have ……………….. value.

Question 5.
The terms having same literal factors are called ………………

Question 6.
The degree of -1x² is ………………… [-2/2]

Question 7.
The degree of 1 is ………………… [1/0]

Question 8.
The degree of the expression 5x – 1 is …………………

Question 9.
The degree of the expression 1 – 3x² is …………………

Question 10.
The degree of the expression 2x² – xy² + x is ………………… [3/6]

Question 11.
The degree of the expression 5 – x + x²y² – 3x³ is ………………… [3/4]

Question 12.
The numerical coefficient of $$\frac {-3}{2}$$ x is ……………….

Answer: $$\frac {-3}{2}$$

Question 13.
The degree of 8 is …………………

Question 14.
In the algebraic expression (3a + 4b – c – 4a + 3ab), the like terms are 3a and ………………..

Question 15.
In the algebraic expression (3a²b + 3ac² + 5ab² – 7a²b² + 5c²a – 20b²a), the like terms are 3ac², 5c²a; 5ab² and ………………….

Question 16.
(2x – 7) is a ………………… [binomial/monomial]

Question 17.
5x is a ……………….. [binomial/monomial]

Question 18.
A term of an expression having no literal factor is called a …………………… term.

Question 19.
An algebraic expression having only one term is called a ………………….

Question 20.
When -3x² is subtracted from 0, we get ……………….

Question 21.
We get ……………….. when (x – 2x²) is subtracted from 1.

Answer: 1 – x + 2x²

Question 22.
The Value of $\frac {x+y}{5}$ at x = 1 and y = 4, is …………………

Question 23.
The value of 4(a + b + c) at a = -1, b = 2 and c = 1, is ……………….