text:
# GE with an intermediate good

intro

I'm looking at a simple model with 1 consumer, 2 goods and 2 firms.

I'm trying to get a price vector [p0, p1] that makes it work.

By makes it work, I mean, supply = demand in all 3 markets.

the problem

The problem is that I'm actually getting a set of price vectors that work.

Consider the picture below :

little details

consumer

There is one consumer who owns both firms and their profits and they go :

$$u = \ln x + \gamma \ln b$$ where $$b$$ is leisure, and $$\gamma$$ is their relative want of leisure

And their budget is :

$$M = wL + \pi_0+ \pi_1$$ where $$w$$ is wage, $$L$$ is their time endowment ($$n + b = L$$), and $$\pi_0$$ and $$\pi_1$$ are each firm's profits.

Solving that gives $$x$$ and $$b$$, so labor supply is $$n = L - b$$.

$$b=\dfrac{\gamma}{1+\gamma} \cdot \dfrac{M}{w}, 0 < b < L$$ $$x=\dfrac{1}{1+\gamma} \cdot \dfrac{M}{p_1}$$ $$n=L-b$$

firms[0]

The first firm, firms[0], uses just labor to make an intermediate good :

$$\pi_0 = p_0 \cdot z_0^{\alpha} - w z_0$$ where $$0 < \alpha < 1$$, $$z_0$$ is their labor demand, and $$y_0 = z_0^{\alpha}$$ is their output.

firms[1]

The other firm, firms[1], uses labor and firms[0]'s output, $$y_0$$.

$$\pi_1 = p_1 \cdot z_1^{\beta} \cdot (k_1+1)^{1-\beta} - w \cdot z_1 - p_0 \cdot k_1$$

Here, $$z_1$$ is their labor demand, $$k_1$$ is their intermediate goods demand, and $$0 < \beta < 1$$. Their output is $$y_1=f_1(z_1,k_1)=z_1^{\beta} \cdot (k_1+1)^{1-\beta}$$

This firm is almost constant returns to scale. That is, if it can make profit at some level $$(z_1,k_1)$$, then it will make more profit at $$(az_1,ak_1)$$ where $$a>1$$. That means that this firm won't settle on some profit maximizing allocation of $$(z_1,k_1)$$ since it will choose to keep buying more, at a given price level. So this firm uses $$y_0$$ as its limiting factor. And from there it decides how much labor to use.

Also, I made it so that $$k_1$$ could be 0. So like, depending on $$[\alpha, \beta, \gamma, L]$$ it could be that firms[0] doesn't even produce.

the markets

So the markets look like this :

Labor : n = z[0] + z[1] @ w

Middle : y[0] = k[1] @ p[0]

Final : y[1] = x @ p[1]

how i'm solving them

First, it's important to note that I do not know what I am doing. I'm doing this purely out of boredom, so please don't be surprised if the answer is something really obvious and I just plum don't know about it.

So I have a little function that takes some random price, say [1, 1], and uses that price vector to get the sum of squares of excess supply, exx. Like, it does "supply minus demand" for each market and squares it, and adds it to the total. It's my way of measuring how bad a price vector is. (Is there a way to measure how bad a price vector is?)

Then it checks a bunch of price vectors around it, like [1+dp, 1], [1-dp, 1], [1, 1+dp] etc.. where dp is the size of the step. And when finds a point around it with a lower exx, it makes that the new price. And repeats. And when it doesn't find a better point, it shrinks dp and does it again.

the problem

The price vector I get changes depending on the starting point. And most of the time I get an exx = 0 (or very very near 0). The problem is that (just based on my graphing it), exx(p[0],p[1]) doesn't seem to be continuous. When I graph exx against p[0] (x-axis) and p[1] (y-axis), I get a whole set (a line) of price vectors and when I check them manually, they work.

centrally planned

When I solve the central planner problem it looks like this :

$$u= \ln x + \gamma \cdot \ln b$$

But $$x=y_1=z_1^{\beta} \cdot (k_1+1)^{1-\beta}$$ and $$k_1=y_0=z_0^{\alpha}$$ so that gives the following :

$$u = \beta \cdot \ln n_1 + (1-\beta) \cdot \ln (n_0^{\alpha} + 1) + \gamma \cdot \ln (L - n_0 - n_1)$$

Take $$\frac{du}{dn_0}$$ and $$\frac{du}{dn_1}$$ and that gives you something like :

$$0 = C_0 \cdot n_0^{1-\alpha} + C_1 \cdot n_0 + C_2$$

Where $$[C_0, C_1, C_2]$$ are constants made from the the exogenous variables $$[\alpha, \beta, \gamma, L]$$

Which is an equation I don't know how to solve but I do it numerically and I can verify that it is the utility-maximizing allocation.

I guess for me the important thing is the solution this gives is on that line I get when I graph $$exx$$ against $$p_0$$ and $$p_1$$. Important insofar as the algorithmic approach I used isn't too wrong, that is.

Anyhow, the reason I want the algorithmic approach to work is because I can easily add lots of consumers, firms, products, firm ownerships etc, and in order to do that, which is fun, I need to make sure my approach, in a technical sense, is sound and actually works.

questions

Does this problem have a unique solution?

Is there a proper way of solving for the solution(s)?

I guess that's all. I can possibly add a link to the work I did.

edit: 2021-06-22

Here's a link to where I'm playing around with this. There seem to be lots of price vectors that work. But since I'm solving them numerically, 'work' really just means 'below a certain level of error'.

https://dactyrafficle.github.io/GE_intermediate/

• Can you write down what the actual production functions are? Also, please use MathJax. Jun 21 '21 at 9:25
• Why do you use $L$ rather than $n$ in the budget constraint? And what does "the markets" part mean? Jun 21 '21 at 19:05
• $L$ goes in the budget constraint because it represents the agent's constraint concerning their leisure/goods allocation. It's not their market income (labor income + dividends), which is $wn+\pi_0+\pi_1$. The markets bit is like this: every input to each market is solely a function of $[w, p_0, p_1]$. So when a given price vector works, supply=demand for each of the 3 markets. I define 'works' as when $\sum (supply-demand)^2$ is 0 or approx 0. Jun 21 '21 at 19:25
• @dactyrafficle I see. So e.g. '@ w' means conditional on the wage? I think you might better clarify it as this way of writing market clear conditions is a little confusing. Jun 21 '21 at 20:11
• @dactyrafficle So you have the FOCs and the market clearing conditions. What prevent you from calculating the equilibrium prices? Jun 21 '21 at 20:24

This is my attempt. The final result is a set of equilibrium equations, which I will not attempt at solving.

### The consumer:

$$\max \ln(x) + \gamma \ln(b) \text{ s.t. } p_1 x + w b = w L + \pi_0 + \pi_1.$$

This is Cobb-Douglass so: \begin{align*} &x_1 = \frac{1}{1 + \gamma} \frac{(wL + \pi_0 + \pi_1)}{p_1},\tag{1}\\ &b = \frac{\gamma}{1 + \gamma} \frac{(wL + \pi_0 + \pi_1)}{w}.\tag{2} \end{align*}

### Firm 0

Firm 1 maxmizes $$\pi_0 = p_0 z_0^\alpha - w z_0.$$ First order condition is: \begin{align*} &p_0 \alpha z_0^{\alpha - 1} = w,\\ \to &z_0 = \left(\frac{w}{\alpha p_0}\right)^{1/(\alpha-1)},\tag{3}\\ \to &y_0 = \left(\frac{w}{\alpha p_0}\right)^{\alpha/(\alpha-1)},\tag{4}\\ \end{align*} Profits are: $$\pi_0 = p_0 \left(\frac{w}{\alpha p_0}\right)^{\alpha/(\alpha-1)} - w\left(\frac{w}{\alpha p_0}\right)^{1/(\alpha-1)}. \tag{5}$$

### Firm 1

The first order conditions for profit maximisation of firm 1 give: $$\beta \frac{y_1}{z_1} = \frac{w}{p_1},\\ (1-\beta) \frac{y_1}{(k_1 + 1)} = \frac{p_0}{p_1}.$$ taking ratios gives: $$\frac{\beta (k_1 + 1)}{(1-\beta)z_1} = \frac{w}{p_0}$$ Now let's substitute this back into the production function: \begin{align*} y_1 &= z_1^\beta (k_1 + 1)^{1-\beta},\\ &= \left(\frac{\beta}{1-\beta}\frac{p_0}{w} (k_1 + 1)\right)^{\beta}(k_1 + 1)^{1-\beta},\\ &= (k_1 + 1) \left(\frac{\beta}{1-\beta}\frac{p_0}{w}\right)^{\beta} \tag{6} \end{align*} Also substituting $$(k_1 + 1)$$ instead gives: $$y_1 = z_1^\beta \left(\frac{1-\beta}{\beta}\frac{w}{p_0}z_1\right)^{1-\beta},\\ = z_1 \left(\frac{1- \beta}{\beta} \frac{w}{p_0}\right)^{1-\beta} \tag{7}$$ Let's use these two expressions to compute the cost function: $$w z_1 + p_0 k_1 = w y_1\left(\frac{1-\beta}{\beta}\frac{w}{p_0}\right)^{1/(1-\beta)} + p_0 y_1 \left(\frac{\beta}{1 - \beta}\frac{p_0}{w}\right)^{1/\beta} - p_0$$ Notice the $$-p_0$$ at the end. So the profit is given by: $$p_1 y - y_1\left(w \left(\frac{1 - \beta}{\beta} \frac{w}{p_0}\right)^{1/(1-\beta} + p_0 \left(\frac{\beta}{1-\beta} \frac{p_0}{w}\right)^{1/\beta}\right) + p_0$$

There is a linear term in $$y_1$$. If the coefficient of this term is strictly positive, then the firm will set $$y_1 = \infty$$. If the term is negative then $$y_1 = 0$$, but this corresponds to an amount $$k_1 = -1$$, which is impossible. So the only choice is where the term in front of $$y_1$$ is zero. This gives: $$p_1 = \left(w \left(\frac{1 - \beta}{\beta} \frac{w}{p_0}\right)^{1/(1-\beta} + p_0 \left(\frac{\beta}{1-\beta} \frac{p_0}{w}\right)^{1/\beta}\right) \tag{8}$$ Then profits of firm 1 are given by $$\pi_1 = p_0 \tag{9}$$

### Equilibrium

On the final goods market, we must have that using $$(1)$$: \begin{align*} &x = y_1,\\ \to&\frac{1}{1 + \gamma} \frac{(wL + \pi_0 + \pi_1)}{p_1} = y_1 \tag{eq-1} \end{align*} On the intermediate goods market, we have using $$(4)$$ and $$(6)$$: \begin{align*} &y_0 = k_1,\\ \to &\left(\frac{w}{\alpha p_0}\right)^{\alpha/(\alpha-1)} = y_1 \left(\frac{\beta}{1-\beta}\frac{p_0}{w}\right)^{1/\beta}-1 \tag{eq-2} \end{align*} On the labour market we have using $$(2)$$, $$(3)$$ and $$(7)$$: \begin{align*} &L - b = z_0 + z_1,\\ \to &L - \frac{\gamma}{1 + \gamma} \frac{(wL + \pi_0 + \pi_1)}{w} = y_1\left(\frac{1- \beta}{\beta} \frac{w}{p_0}\right)^{1/(1-\beta)} + \left(\frac{w}{\alpha p_0}\right)^{1/(\alpha-1)}, \tag{eq-3}\\ \end{align*} We also have, from $$(5)$$ and $$(9)$$ the profits of both firms: \begin{align*} &\pi_1 = p_0, \tag{eq-4}\\ &\pi_0 =p_0 \left(\frac{w}{\alpha p_0}\right)^{\alpha/(\alpha-1)} - w\left(\frac{w}{\alpha p_0}\right)^{1/(\alpha-1)} \tag{eq-5} \end{align*} and finally, the equilibrium condition for firm 1: $$p_1 = \left(w \left(\frac{1 - \beta}{\beta} \frac{w}{p_0}\right)^{1/(1-\beta} + p_0 \left(\frac{\beta}{1-\beta} \frac{p_0}{w}\right)^{1/\beta}\right) \tag{eq-6}$$ Let's count unknowns:

• prices $$w, p_0, p_1$$

• output level $$y_1$$

• profits $$\pi_1, \pi_0$$.

we have 3 market equilibrium conditions, 2 profit equations and 1 pricing equation for firm 1. So in total 6 equations in 6 unknowns, which should allow us to solve the system.

Note: Not all equations are independent as equilibrium on two markets will automatically give equilibrium on the third one. However, we can normalize price of one good to unity, so we end up with 5 equations and 5 unknowns.

• The way you wrote the consumer's budget constraint is good. In $(3)$ I believe you're missing an $\alpha$ in the denominator. $(6)$ to $(9)$ is insightful. My thinking was that the firm would want to buy as much $z$ and $k$ along the expansion path, and when it ran out of one, it would take the limiting input as fixed, and then do a 1 variable max problem. $(9)$ is very interesting too, because when I solve the central planner problem numerically, and reverse-engineer it to get prices that make the market work at that allocation, that is the result I get too. Jun 22 '21 at 19:23
• @dactyrafficle I fixed the mistake (I hope).
– tdm
Jun 23 '21 at 5:24
• I'm a little bit intrigued to the logic of $(6)$ thru $(9)$. A minor correction: If you write $A=\frac{(1-\beta) \cdot w}{\beta \cdot p_0}$ then $(8)$ should be $wA^{\beta-1}+p_0A^{\beta}$. Solving that let's me complete the price set $[w, p_0,p_1]$ that coincides with the centrally planned allocation. In fact, the solution is part of the set of solutions I get when I use a numerical approach. I'm still looking at it but this so far I think it pretty insightful. Jun 24 '21 at 3:08


text:
# Is “sensitivity” a term in numerical computation?

The section 4.2 "Poor Conditioning" in the book Deep Learning defines the condition number of the function $$f(x) = A^{-1}x$$ as

\begin{align} \underset{i,j}{\max}~ \Bigg| \frac{\lambda_i}{ \lambda_j} \Bigg|. \end{align}

which is the ratio of the magnitude of the largest and smallest eigenvalue.

and says

When this number is large, matrix inversion is particularly sensitive to error in the input. This sensitivity is an intrinsic property of the matrix itself, not the result of rounding error during matrix inversion.

I guess "sensitivity" here means how rapidly a function changes with respect to small changes in its inputs.

Is "sensitivity" a term, something like the one in Medical statistics or plain English here?

Yes, its commonly used to describe the magnitude of an outputs change with respect to an input. In CFD optimization we provide a sensitivity vector to the optimizer (or the derivatives of each output to each input). In this case the matrix inversion is sensitive to (for example) small differences on the right hand side or floating point error.

• thank you. Does "CFD" here refers to "Computational fluid dynamics"? – shi95 May 1 at 3:33
• @EMP Is there some difference with the "conditioning" of a problem? From what you described I immediately think about well/ill conditioned problems – VoB May 1 at 8:17
• @shi95, yes it does. – EMP May 1 at 19:41
• @VoB, I'm sorry I don't exactly understand your question. – EMP May 1 at 19:41
• @EMP I think Federico's post is the answer to my comment – VoB May 2 at 16:17

It is "a term" that is used commonly (it appears many time, for instance, in Higham's classical book Accuracy and Stability of Numerical Algorithms), but it is used with its plain English meaning, in my experience. You can use it in a sentence like

The condition number measures the sensitivity of $$A^{-1}$$ with respect to changes in $$A$$.

but it does not have a precise mathematical definition like "condition number" or the one in statistics.

format: plain


text:
# I need help converting Unity .shader to Shader Graph: Bumpy Glass

I've been tasked with converting an older Unity shader which does Bumpy Glass with some alpha modification, to URP. I'm assuming that means using Shader Graph. I'm a bit new to shaders in general, but managed to find various clues here and there and have come to the conclusion that creating a Custom Function is probably the easiest way (I could be wrong). Did I mention, this is targeted for WebGL?

Shader "Custom/Glass-Seethrough Alpha" {
Properties{
_Color("Main Color", Color) = (1,1,1,1)
_SpecColor("Specular Color", Color) = (0.5, 0.5, 0.5, 1)
_Shininess("Shininess", Range(0.03, 1)) = 0.3
_MainTex("Base (RGB)", 2D) = "white" {}
_BumpMap("Normalmap", 2D) = "bump" {}
_AlphaMap("Alphamap", 2D) = "black" {}
_DistAmt("Distortion", range(0,256)) = 100
}
GrabPass { }

Tags { "Queue" = "Transparent" "IgnoreProjector" = "True" "RenderType" = "Opaque" }
LOD 200

CGPROGRAM
#pragma target 3.0
#pragma exclude_renderers gles
#pragma vertex vert
#pragma surface surf BlinnPhong alpha
#include "UnityCG.cginc"

float4 _Color;
float _Shininess;
sampler2D _MainTex;
sampler2D _BumpMap;
sampler2D _AlphaMap;
sampler2D _GrabTexture;
float _DistAmt;
float4 _GrabTexture_TexelSize;

struct Input {
float2 uv_MainTex;
float2 uv_BumpMap;
float2 uv_AlphaMap;
float4 proj : TEXCOORD;
};

void vert(inout appdata_full v, out Input o) {
UNITY_INITIALIZE_OUTPUT(Input, o);
float4 oPos = UnityObjectToClipPos(v.vertex);
#if UNITY_UV_STARTS_AT_TOP
float scale = -1.0;
#else
float scale = 1.0;
#endif
o.proj.xy = (float2(oPos.x, oPos.y * scale) + oPos.w) * 0.5;
o.proj.zw = oPos.zw;
}

void surf(Input IN, inout SurfaceOutput o) {
half3 tint = tex2D(_MainTex, IN.uv_MainTex);
half3 nor = UnpackNormal(tex2D(_BumpMap, IN.uv_BumpMap));
half3 alpha = tex2D(_AlphaMap, IN.uv_AlphaMap);

float2 offset = nor * _DistAmt * _GrabTexture_TexelSize.xy;
IN.proj.xy = offset * IN.proj.z + IN.proj.xy;
half4 col = tex2Dproj(_GrabTexture, UNITY_PROJ_COORD(IN.proj));
half3 tinted = col * tint * _Color;

o.Albedo = lerp(col, tinted, alpha);//tex.rgb * _Color.rgb * col.rgb;
o.Normal = nor.rgb;
o.Specular = _Shininess;
o.Gloss = col.a;
o.Alpha = alpha.r;
}
ENDCG
}
FallBack "Diffuse"

}

The idea is to have the bump map animate on, say, uv.v, to simulate water flow through a pipe, and animate the alpha map to 'reveal' the pipe. Here's what it should look like:

I found this article which helped (I THINK!) with the 'vert' part of the shader. Where I'm getting stuck is in the 'surf' part. I've set up a bunch of inputs, and maybe I've gone down the wrong rabbit hole, but I'm getting an error saying that uv_MainTex does not exist:

Here's a shot of the work in progress:

And the full text (admittedly incomplete):

half3 tint = SAMPLE_TEXTURE2D(MainTex, SS, uv_MainTex);
half3 nor = UnpackNormal(tex2D(BumpMap, uv_BumpMap));
half3 alpha = tex2D(AlphaMap, uv_AlphaMap);

float2 offset = BumpMap * DistAmt; * GrabTexture_TexelSize.xy;
proj.xy = offset * proj.z + proj.xy;
half4 col = tex2Dproj(GrabTexture, proj);
half3 tinted = col * tint * color;

Albedo = lerp(col, tinted, alpha);

But I'm getting this error that uv_MainTex does not exist.

Am I even on the right track? Am I using the wrong macros? Do I need additional SamplerStates? Or hey, has someone already made such a shader I can get and call it a day? Thanks

New contributor
Dave is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.

This should give you something fairly close to what you had. The proj data coming from the vertex shader wasn't doing anything geometrically meaningful, so I skipped that for simplicity. You can get a similar effect by tuning your Distortion parameter instead.


text:
# Learning Userland DRM

## Part I

### January 23, '21

DRM is the bridge between the Linux kernel and actual software running atop Linux, for most operations which involve a graphics card. For example, X11 and Wayland actually send drawing commands to the graphics card using DRM.

In this post I’ll be talking about how to write software which communicates with the graphics card directly, completely bypassing X. Do mind that I am also learning while I write, so some information might be inaccurate and also I’m jumping through unnecessary hoops, from the perspective of just getting things done.

I’ll be using C++.

# libdrm

The kernel DRM API is file based, which means you can open a file representing the graphics card (usually /dev/dri/card0) and use that file to send actual commands to the card, using the ioctl() syscall with the opened file representing the card.

Fortunately, there is also a library conveniently named libdrm that simplifies all that. Keep in mind that all that it does is what we mentioned above: opening a file and sending commands through ioctl() syscalls.

In C++ a simple program that interacts with the video card through DRM might look like this at the top level:

#include <cassert>
#include <unistd.h>
#include <fcntl.h>

// The file which represents the card
static const char gpu_dev[] = "/dev/dri/card0";

int open_gpu() {
return open(gpu_dev, O_RDWR);
}

int close_gpu(int gpu) {
return close(gpu);
}

int main(int arg_count, char** arg_vector) {
const int gpu = open_gpu();
assert(gpu >= 0);

// some ioctl() magic here...

assert(close_gpu(gpu) >= 0);

return 0;
}

And as a sanity check we can try to run it in a Linux machine:

❯ clang++ drm-part-1.cpp -o drm-part-1
❯ ./drm-part-1
❯

Surely enough, it runs! But it doesn’t do anything yet. It needs that ioctl() magic, but at least no assertions failed and we know we can open the video card file now!

Lets include another library which has some very handy structs and also definitions of parameters for ioctl(): libdrm/drm.h. While we are at it, lets also include sys/ioctl.h and iostream for syscalls and printing.

#include <libdrm/drm.h>
#include <sys/ioctl.h>
#include <iostream>

# Driver Version

The first interaction that we are going to have with the card is not going to be with the card itself but just with the driver (bummer, I know). Lets try to retrieve the driver version using a ioctl() call.

There is a definition for a struct which holds exactly this data in libdrm/drm.h, drm_version_t:

struct drm_version {
int version_major;	  /**< Major version */
int version_minor;	  /**< Minor version */
int version_patchlevel;	  /**< Patch level */
__kernel_size_t name_len;	  /**< Length of name buffer */
char *name;	  /**< Name of driver */
__kernel_size_t date_len;	  /**< Length of date buffer */
char *date;	  /**< User-space buffer to hold date */
__kernel_size_t desc_len;	  /**< Length of desc buffer */
char *desc;	  /**< User-space buffer to hold desc */
};

To get the version from the driver we need to do three things:

1. Initialize a drm_version_t, notice there are some char * fields for which we need to allocate memory.
2. Call ioctl() with the appropriate parameters and a pointer to our drm_version_t.
3. Print out the information we get and free everything.

Seems easy enough right? Lets do it!

drm_version_t create_version() {
drm_version_t version;

version.name = new char[256];
version.date = new char[256];
version.desc = new char[256];

return version;
}

And:

void delete_version(const drm_version_t &version) {
delete version.name;
delete version.date;
delete version.desc;
}

To print it out:

void print_version(const drm_version_t &version) {
std::cout << std::string(version.name, version.name_len) << std::endl
<< std::string(version.desc, version.desc_len) << std::endl
<< version.version_major << "."
<< version.version_minor << "."
<< version.version_patchlevel << std::endl;
}

And finally, we can call all that in our main(), notice the ioctl() call using the DRM_IOCTL_VERSION value, which is defined in libdrm so we don’t have to worry about the specific values passed:

int main(int arg_count, char** arg_vector) {
const int gpu = open_gpu();
assert(gpu >= 0);

drm_version_t version = create_version();
assert(ioctl(gpu, DRM_IOCTL_VERSION, &version) == 0);
print_version(version);
delete_version(version);

assert(close_gpu(gpu) >= 0);

return 0;
}

After a little compiling, behold the power!

❯ ./drm-part-1
i915
Intel Graphics
1.6.0

We are talking to the DRM driver of my Intel graphics card! No X11 involved, this runs even directly in the tty!

Keep in mind that libdrm does a lot more than just provide the definitions and we do not need to be calling ioctl() directly, but for the purposes of learning I found important know exactly how the syscalls are made to the graphics driver.

drm-part-1.cpp


text:
## Introduction

The Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES) is a binding and multilateral environmental agreement with the objective of prohibiting trade of wildlife species in immediate danger of extinction (Appendix I listing) and preventing additional threatened species reaching this critical stage by regulating trade so that it is traceable, legal, and sustainable (Appendix II listing). CITES is the first line of defense against illegal wildlife trafficking, with 183 parties (i.e. signatory nations), and lists over 35,000 animal and plant species1. The vast majority of these species (96%) are listed in Appendix II and exporting parties have the obligation to document that traded specimens are traceable through the supply chain, were legally obtained, and that trade it is not detrimental for the survival of the species2. Importing parties are required to monitor imports and certify that incoming specimens are accompanied by the appropriate documentation2. Failure to properly implement CITES for any listed species can eventually lead to international trade sanctions being placed on offending parties2.

Sharks are a group of marine fish many of which are threatened by overexploitation to satisfy demand for internationally traded products, primarily fins (for use in luxury Asian soup dishes) and meat3,4. Twelve shark species have been listed in CITES Appendix II to date. The first round of shark listings took place from 2001–2004, with the whale shark Rhincodon typus (2001), basking shark Cetorhinus maximus (2001), and great white shark Carcharodon carcharias (2004). These species all share the characteristics of being iconic, extremely large bodied (>4.0 m total length at maturity), and some products from these species (e.g., fins, dressed carcasses) can be readily identified by their very large sizes (except when procured from juvenile stages) or other morphological features (e.g., jaws and teeth for great whites). These species also share the characteristic that Appendix II listing is not as stringent as the domestic protection these species receive in many jurisdictions (i.e., landing and trade prohibition4,5). Products of these species have generally not been detected in post-listing studies of major shark fin and meat markets, suggesting they are rare in trade6,7,8,9,10.


text:
# Introduction to Complex Systems

Prof. Dirk Brockmann, Winter Term 2021

# Two-dimensional Dynamical Systems

In the last few lectures we discussed one dimensional dynamical systems like this one: [ \dot{x}=f(x) ] and learned how to analyze their asymptotics graphically and by calculation. Essentially, this amounted to just plotting the function $f(x)$ and finding the zeros, see One Dimensional Dynamical Systems.

Some 1d-systems were motivated by underlying reaction kinetic cartoons, for example [ X\xrightarrow{\alpha}2X\qquad X+X\xrightarrow{\beta}\emptyset. ] In this particular example, replication (left reaction) and annihilation (right equation) compete'' and balance. This example lead to the dynamical system [ \dot{x}=\alpha x-\beta x^{2}=x(\alpha-\beta x) ] for the concentration or abundance of $X$-particles.

We discussed that in systems like this one we can schematically understand the system by a drawing like the one in Figure 1 in which each link indicates the type of influences that play a role. In this case reaction (1) is a positive autocatalytic influence of $X$ onto itself and reaction (2) is an inhibitory influence. Just by looking at the diagram we see that these forces'' can balance yielding a stable equilibrium.

If we look at a different dynamical system, for instance [ 3X\xrightarrow{\alpha}4X\qquad X\xrightarrow{\beta}\emptyset ] with a corresponding dynamical system [ x=x\left(\alpha x^{2}-\beta\right) ] we see that it is structurally is different. However, schematically, it looks the same, we have a reaction that is a positive influence (1) and increases the abundance, and another equation (2) that decreases it. And also here we can expect a balancing scenario in which both reactions stabilize the system.

Sometimes, we can already see what to expect from a dynamical system just by drawing the schematics. This is even more powerful in 2D where analytic treatments are often much harder and developing an intuition for the dynamics even more vital.

## Basics

Now we will discuss two dimensional dynamical systems. Don't worry, we won't be increasing dimensionality step by step to discuss peculiarities of each of the dimensions. We will only discuss two dimensions in more detail because even higher dimensional systems can sometimes be understood by putting 2D-motifs together that are derived from 2D-interactions. We'll discuss this in more detail later.

A two dimensional system's state is defined by two dynamic variables, usually denoted by $x$ and $y$ so the state is the pair: [ \mathbf{r}(t)=\left(x(t),y(t)\right). ] As you may expect, we will be discussing systems that are defined by [ \dot{x} =f(x,y) ] [ \dot{y} =g(x,y), ] so by two functions $f$ and $g$ that define how each component of the state changes. A more compact notation is this one [ \mathbf{\dot r}=\mathbf{f}(\mathbf{r}) ] where the $\mathbf{r}(t)$ is a 2D state-vector and a state-dependent vector [ \mathbf{f}(\mathbf{r})=\left(f(x,y),g(x,y)\right) ] that tells us how the state changes in every state $\mathbf{r}$. This notation becomes more practical in higher dimensions. For 2D we can stick with the explicit notation for each component using $x$ and $y$. If we visualize the state-space as the $x-y-$plane, then we can think of $\mathbf{f}(\mathbf{r})$ as little vectors defined at every location which is why this is called a vector-field. We'll discuss this in more detail below.

## An example

First let's derive a 2D dynamical system using reactions and schematics. Let's say we have two different types of particles, $A$ and $B$ and they interact with each other and with themselves according to four reactions like so:

[ A\xrightarrow{\alpha}2A \qquad A+A\xrightarrow{\beta}A ]

[ B\xrightarrow{\gamma}\emptyset\qquad A+B\xrightarrow{\delta}2B ]

Let's denote the concentration of $A$ and $B$ by $x$ and $y$, respectively. Then we know that reaction (1) above replicates $A$ and therefore increases $x$ by a term $\alpha x$. The second equation models a collision of $A$'s by which one agent is destroyed, so this is a term $-\beta x^2$. The third reaction decrease $B$ by spontaneous disappearance so we have a terms $-\gamma y$ by which $y$ change accordingly.

The interesting equation is (4) which is a two particle reaction of $A$ with $B$ that both decreases $x$ and increases $y$ by a factor $\delta xy$. So, all in all we get a dynamical system

[ \dot{x} =\alpha x-\delta xy-\beta x^2 ] [ \dot{y} =\delta xy-\gamma y ] Schematically we see that $A$ is self-reinforcing, both $A$ and $B$ have an inhibitory impact on themselves (reactions (2) and (3)). The interesting part is the interaction that is both, a positive impact of $A$ onto $B$ and a negative impact of $B$ onto $A$ which is driven by the last reaction. So, the interaction is of opposite sign in a way. Figure 2 illustrates the schematic diagram for this system.

## Schematic analysis

In the example above we didn't do much, except made a little schematic. We can do more. We can compute fixpoints, their stability and gain an overall understanding of the behavior and of the asympotics of a system. Below is a recipe on how to investigate a system like the example system above. Before we do this for some examples let's see however what kind of qualitatively different schematics we can have in a 2D-system.

### Mutual Inhibition

Let's look at these reactions

[ \emptyset \xrightarrow{\alpha}A \qquad A \xrightarrow{\beta}\emptyset ]

[ B \xrightarrow{\gamma} 2 B \qquad 3B \xrightarrow{\delta}2B ]

[ A+B\xrightarrow{\xi}\emptyset ]

The first two equations alone yield a dynamical systems [ \dot{x}=\alpha-\beta x ] so just these two would imply a self-regulating $A$ population by spontaneous creation of $A$'s (reaction 1) and spontaneous death (reaction 2). So by itself we have something that calls for an equilibrium in which both reactions balance. In fact, this is what you get if you just have $A$'s. If we only focus on $B$ reactions (3) and (4) we get [ \dot{y}=\gamma y-\delta y^{3}. ] Again a system in which two forces'' of opposite sign balance. The interesting stuff is happening because of equation (5) which implies a negative change of both $A$ and $B$ proportional to $-\xi xy.$ This couples the system and we get

[ \dot{x} =\alpha-\beta x-\xi xy ] [ \dot{y} =\gamma y-\delta y^{3}-\xi xy ] Unlike the system above, schematically we see a system in which the interaction is mutually inhibitory.

So by looking at the corresponding schematic in Figure 3 we can imagine that in such a scenario we could have an asymptotic state in which either $A$ and $B$ win'' and push the competitor to zero, or we may expect that for a low value of $\xi$ in a regime of weak inhibition that both $A$ and $B$ approach a non-zero stable state each.

You can investigate and confirm this conjecture using the Interactive tool for 2D systems in the Script Collection. In the tool, choose the Mutual Inhibition system and tune the parameters to generate a stable coexistence state for some set of parameters, and try to find conditions for a winner takes it all scenario.

### Cooperative systems

Sometimes, we obtain the ODEs defining a dynamical system or a model in a different way, e.g. a book, or a paper and it helps reading a schematic interaction off the equations defining the system. Consider this system [ \dot{x} =r_{1}x(1-x)-x+\alpha xy(1-y) ] [ \dot{y} =r_{2}y(1-y)-y+\beta yx(1-x) ] defined for $x,y\geq\text{0}$ and say for parameters $0\leq r_{1/2}<1$. Let's clarify this equation a bit and write it as:

[ \dot{x} =r_{1}\,\lambda(x)-x+\alpha x\,\lambda(y) ] [ \dot{x} =r_{2}\lambda(y)-y+\beta y\,\lambda(x) ] where the function [ \lambda(z)=z(1-z) ] is the typical, self-regulated logistic growth function that is linear for low concentrations and has the limiting $1-z$ factor. We see that the interaction is hidden in the last terms of the rhs of the two parts of the system and controlled by the two parameters $\alpha$ and $\beta$. Without the interaction terms each of the systems is equivalent to the epidemic SIS system with a subcritical reproduction number (see Script on One-Dimensional Dynamical Sysems). So for instance if $\alpha=0$ the first equation is [ \dot{x}=r_{1}\,x(1-x)-x ] that only has the $x=0$ state as a stable fixpoint (because we say that $r_{1}<1$) and likewise for $y$ if $\beta=0$. Now let's look at what the term $\alpha xy(1-y)$ does. Well this term is positive so it's a cooperative influence of $y$ onto $x$. But it's limited if $y$ gets too large the factor $1-y$ makes sure the cooperativity is clamped. The same is true the other way around. The term $\beta yx(1-x)$ represents a positive effect of $x$ onto $y$. So if we were to draw a schematic we have a system in which each element is doing the usual self-excitation with saturation, and additional negative self-loop but a positive interaction like so:

A system like this is mutually cooperative. The questions now are: What is the asymptotic behavior of the system? Is the cooperativity strong enough to stabilize a stationary state in which both $x$ and $y$ are non-zero? Is the $x=y=0$ state always stable? Do we have states in which only one species is stabilized by the other, but not vice versa?

You can see for yourself using the Interactive tool for 2D systems again. In the systems menu pick Cooperative System. It turns out it can can up to 5 different stationary states and three stable coexistence states. Try to find conditions for it.

### Activator-inhibitor systems

The first example we discussed and schematically analyzed is an example of an activator-inhibitor system in which one element inhibits the other and the other activates the first. These systems are defined by an asymmetric coupling. One of the more famous systems that captures activator-inhibitor dynamics is a predator-prey model known as the....

#### ...Lotka Volterra model

This is defined by the equations

[ \dot{x} =x(\alpha-\beta y) ] [ \dot{y} =y(\gamma x-\delta) ] for the abundance $x$ and $y$ of two species $A$ and $B$ respectively. Species $B$ inhibits $A$, see the term $-\beta xy$. The reproduction of $B$ increased by the abundance of $A$, see the term $\gamma xy$. Therefore, $A$ must be the prey that get consumed by the predator $B$. In addition to the interaction terms we have regular replication in the prey (the term $\alpha x$) and decay in the predator $-\delta y$. This therefore follows a similar schematic as shown in Figure 2. Note, however, that in the absence of the predator, prey will grow exponentially without bounds, so the LV-system makes no sense in this marginal case. Something to keep in mind.

In an activator-inhibitor system like this one, we can imaging that the interaction of activation and inhibition can generate a stable state in which both effects balance. But we can also imagine that the system cycles through phases:

1. If both predator ($y$) and prey ($x$) are low, $x$ can exponentially increase.
2. After that, the predator abundance can pick up because at some point there's sufficient prey to consume.
3. The increase in predators will decrease in prey until ...
4. .... there's not enough prey and the predator population will decrease again because there's not enough food, yielding situation 1.

We see just by this train of thoughts that we can expect activator-inhibitor systems to exhibit cyclic solutions.

## Analysis

One way of investigating 2D systems is of course just by plotting $x(t)$ and $y(t)$ for an initial condition $x(0)$ and $y(0)$ as a function of $t$. But just like in 1D systems, if we are more interested in the asymptotics, it's more helpful to visualize the trajectory in phasephase, sometimes called state-space, which is the 2D $x-y$-plane. The initial condition is just a location in this space, and the solution ($x(t),y(t)$) is a path in it. Irrespective of what this solution looks like we can deduce a few things that must hold. For example, trajectories can't cross. If this were possible, the solution would not be unique at the intersection. It wouldn't know where to go.

## The flow field

As mentioned above we should think of a 2D system as a vector-field, so at every position ($x,y$) in state-space we imagine a vector [ \left(f(x,y),g(x,y)\right) ] Because in the dynamical system the components of this vector represent the change of $x$ and $y$ the vectors point in the direction of motion. Now we can put a grid across a region of state-space and at each grid point draw a vector. This is difficult to do by hand, but computers can do this very quickly. See Panel XX. In the panel the, the flow field of the cooperative system

[ \dot{x} =r_{1}\,x(1-x)-x+\alpha x\,y(1-y) ] [ \dot{y} =r_{2}y(1-y)-y+\beta y\,x(1-x) ] is shown for $r_{1}=0.1$, $r_{2}=0.1$, $\alpha=3.974$ and $\beta=4.641$ .

Just by looking at the little vectors we see that potentially the system has 3 attractors which are points in state-space. In the panel you can move the mouse across the state-space and for an initial condition under the mouse pointer'' a full trajectory is computed. When you click you can freeze it and select another one. Try it.

## Asymptotics Brainstorm

We are almost ready to do some math. But before we do, let's take the previous example and discuss what we could expect from a 2D dynamical system when we think of asymptotics.

Well, we could have points as attractors, fixpoints in 2D. So points like $(x^{\star},y^{\star})$ that not only fulfill [ f(x^{\star},y^{\star})=g(x^{\star},y^{\star})=0 ] but that also attract trajectories. In the example above we observe three of those attractive fixpoints and depending on the initial condition different trajectories will approach one or the other fixpoint.

We can also think of other types of attractors in 2D. In two dimensions we can imaging that trajectories can be closed loops, cyclic solutions such that [ x(t) =x(t+T)\qquad y(t) =y(t+T) ] for any point on the cycle and $T$ being the time it takes to go around. This is the sort of thing we expect for activator-inhibitor systems. And in fact we will discuss and analyze a specific example of this. You can sneak-peak in the Examples of 2D Dynamical Systems - Interactively - Script and see that the Lotka-Volterra system actually exhibits a family of such periodic orbits around one of its fixpoints.

### Poincare-Bendixon Theorem

In addition to this we could have trajectories that explode to $\pm\infty$. Let's exclude this for a second and assume that we can define a region in state-space for which we can show that the flow enters the region everywhere on its boundary. What can we say about trajectories that enter that region. For sure there could be an attractive fixpoint in the region that attracts all trajectories that enter the region. We could also have a limit cycle, a periodic solution that attracts all the trajectories, see for example the Brusselator or Predator Prey system in the Examples of 2D Dynamical Systems - Interactively - Script . But could we also have something like a trajectory that folds up infinitely many times like illustrated in Figure 5?

Could we have a trajectory that folds itself ad-infinitum, never crossing but getting folded closer and closer? In principle this should be possible. But it turns out that a deep mathematical theorem known as the Poincare-Bendixon theorem states that in 2D this isn't possible. So all we are left with as asymptotic behaviors are periodic cycles or stable fixpoints (potentially infinitely many, though). This limits the variety of possibilities immensely. The Poincare-Bendixon Theorem is only true for 2D. When we have an additional dimension, we can have bounded solutions that are neither periodic not fixed. We can have what is known as strange attractors. These objects are very weird and play a role in deterministic chaos.

## Fixpoints and nullclines

Now let's dig into 2D systems a bit deeper. Given a system we usually follow a recipe. First, just like in 1D we try to find the systems fixpoints. These are all the points in state-space for which

[ \dot{x} =f(x,y)=0 ] [ \dot{y} =g(x,y)=0 ] Let's do this for the Lotka-Volterra system: [ \dot{x} =x(\alpha-\beta y) ]

[ \dot{y} =y(\gamma x-\delta) ]

This yields the trivial fixpoint [ x^{\star}=y^{\star}=0 ] and the second fixpoint [ (x^{\star},y^{\star})=(\delta/\gamma,\alpha/\beta) ] First thing we do is draw those in the state space. For the Lotka-Volterra system we are only interest in the upper-right quadrant of the $x-y-$plane, i.e. the region for which $x\geq0$ and $y\geq0$.

The next thing we do is draw the nullclines. The $x$-nullcline is defined by all the points that fulfill [ f(x,y)=0 ] Usually this is a set of simple curves defined in state space. For the LV-system the $x-$nullclines are the points defined by $x=0$ which is the $y$-axis and the horizontal line [ y=\alpha/\beta ] Because $\dot{x}=0$ on these curves the trajectories of the LV system must cross the $x-$nullcline vertically, i.e. parallel to the $y$-axis. The $x$-nullclines also separate the state-space into regions where $\dot{x}>0$ (below the nullcline $\alpha/\beta$) and where $\dot{x}<0$, above that line. Likewise the $y$-nullclines are given by the $x-$axis (where $y=0$) and the vertical line $x=\delta/\gamma$. The trajectories of the system have to cross the $y$-nullcline horizontally. The $y$-nullcline separates the state-space into regions where $\dot{y}>0$ and $\dot{y}<0$.

Fixpoints are obviously the intersections of $x$- and $y$-nullclines. For the Lotka-Volterra system we can separate the state-spaces into four region in which the sign of the rate of change $\dot{x}$ and $\dot{y}$ are different. So here we have 4 regions:

• A: $\dot{x}<0$ and $\dot{y}<0$ so trajectories point south-west
• B: $\dot{x}>0$ and $\dot{y}<0$ so trajectories point south-east
• C: $\dot{x}>0$ and $\dot{y}>0$ so trajectories point north-east
• D: $\dot{x}<0$ and $\dot{y}>0$ so trajectories point north-west

So by drawing little arrows in these directions we can get a sense of the trajectories. It looks like that, as expected for an activator-inhibitor system that the trajectories circle around the fixpoint. This is also clear from looking at the flow-field.

All of this is information is collected in Figure 6.

One of the things we cannot say here whether the circular movement in state-space will spiral down to the fixpoint in the center, or spiral outwards or generate concentric periodic orbits or potentially have a limit cycle. It turns out that for the LV system, we can show that all the orbits are periodic, the shape depending on the initial condition. You can check this again here.

## Fixpoint stability

To get a bit more information about a system we can study the behavior of a system near its fixpoints. Remember that in 1D $\dot{x}=f(x)$ with a fixpoint at $x^{\star}$ it was sufficient to look at the derivative at the fixpoint $df(x)/dx\,|\,x=x^{\star}$. Likewise in a 2D system, we can linearize the system around the fixpoint, Let's say [ x(t) =x^{\star}+\delta x(t) ] [ y(t) =y^{\star}+\delta y(t) ]

where $\delta x,\delta y$ terms are small deviations from the fixpoint. So this yields [ \dot{\delta x} =\dot{x}=f(x^{\star}+\delta x(t),y^{\star}+\delta y(t)) ] [ \dot{\delta y} =\dot{y}=g(x^{\star}+\delta x(t),y^{\star}+\delta y(t)) ] which is approximately [ \dot{\delta x} =\dot{x}\approx f(x^{\star},y^{\star})+ \frac{\partial f(x^\star,y^\star)}{\partial x}\delta x + \frac{\partial f(x^\star,y^\star)}{\partial y}\delta y ] [ \dot{\delta y} =\dot{y}\approx g(x^{\star},y^{\star})+ \frac{\partial g(x^\star,y^\star)}{\partial x}\delta x + \frac{\partial g(x^\star,y^\star)}{\partial y}\delta y ] The first terms in each equation are zero. So we can write this compactly as

$$\dot{\left(\begin{array}{c} \delta x\\ \delta y \end{array}\right)}=\mathbf{J}(x^{\star},y^{\star})\left(\begin{array}{c} \delta x\\ \delta y \end{array}\right)$$

where the matrix $\mathbf{J}$ is the Jacobian, the matrix of partial derivatives at the fixpoint

$$\mathbf{J}(x^{\star},y^{\star})=\left(\begin{array}{cc} \partial_{x}f & \partial_{y}f\\ \partial_{x}g & \partial_{y}g \end{array}\right)_{(x,y)=(x^{\star},y^{\star})}$$

Now it turns out that the stability of the fixpoint is determined by the eigenvalues of the Jacobian. This is why. The dynamics for the variables $\delta x$ and $\delta y$ is linear and has the form

[ \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} ]

if we write the two component vector compactly as $\mathbf{x}$ and $\mathbf{A}$ is a $2\times 2$ matrix

$$\mathbf{A}=\left(\begin{array}{cc} a & b\\ c & d \end{array}\right)$$

A solution to the above system is (generically) given by

[ \mathbf{x}(t)=c_{1}\mathbf{v}_{1}e^{\lambda_{1}t}+c_{2}\mathbf{v}_{2}e^{\lambda_{2}t} ]

where $c_{1/2}$ are constants, $\mathbf{v}_{1/2}$ the eigenvectors of $\mathbf{A}$ and $\lambda_{1/2}$ the eigenvalues.

Also, if your imagine that $\mathbf{A}$ has two different and linearly independent eigenvectors $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$ pointing in some direction and if you start your system on the line that is spanned by an eigenvector the trajectory will stay on it.

So along that line the system behaves like a one-dimensional system. And we can have different combinations of stable or unstable lines that are spanned by the two eigenvalues. So just by this reasoning we can expect three different scenarios, see Figures 7-9.

The situation is a bit complexified by the fact that the eigenvalues and eigenvectors of a general $2\times 2$ matrix can be complex numbers like $z=\alpha+i\beta$ where $i^2=-1$ is the imaginary unit. If you investigate linear systems comprehensively you can dig into this and classify all possible types of behaviors of linear systems. What matters here is that we can see that the asymptotic behavior is generically dominated by the real part of the eigenvalues. We can see already that if one of the real parts is positive one of the terms will explode to infinity which means that the fixpoint is unstable. Stability requires that the real part of both eigenvalues are non-positive.

Without going into the dynamics of linear systems too deeply, let's briefly summarize what can happen. To this end we need to compute the eigenvalues of a general $2\times 2$ matrix. These can be expressed in terms of the trace $s$ and the determinant $\Delta$ of the matrix which are defined as [ s =a+d ] [ \Delta =ad-bc ] The eigenvalues are then [ \lambda_{1/2}=s/2\pm \left(s^{2}/4-\Delta \right)^{1/2} ] It's good practice to check this. We can now classify different situation.

### Saddle nodes: $\Delta<0$

In this case both eigenvalues are real, one is negative one is positive. This means the fixpoint is a saddle with one direction that is attractive and one direction that is repelling. Trajectories around such a fixpoint look like the ones in .

### Stable node: $s^{2}/4>\Delta>0$ and $s<0$

In this case both eigenvalues are real and smaller than zero. So we have two directions that are both attractive. and all initial conditions in the vicinity of the fixpoint will move towards it.

### Unstable node: $s^{2}/4>\Delta>0$ and $s>0$

In this case both eigenvalues are real and larger than zero. So we have two directions that are both repulsive. All initial conditions in the vicinity of the fixpoint will move away from it.

### Stable spiral: $\Delta>s^{2}/4$ and $s<0$

Because of the first condition, the eigenvalues are complex and also complex conjugates of each other. [ \lambda_{1/2}=s/2\pm i(\Delta-s^{2}/4)^{1/2} ]

This implies an oscillatory component of the solution. The term under the root is the angular frequency with with a solution oscillates around the fixpoint. Because of $s<0$ the trajectories approach the fixpoint and it is therefore a stable spiral.

### Unstable spiral: $\Delta>s^{2}/4$ and $s>0$

This is just the opposite scenario. We still have a spiral but it expands outwards and is an unstable spiral.

### Focus: $s=0$ and $\Delta>s^{2}/4$

In this case the eigenvalues are purely imaginary and trajectories of the linear system are periodic orbits that cycle around the fixpoint.

This is a non-comprehensive list of behaviors, there are other special cases for instance when $\Delta=s^{2}/4$ or $\Delta=0$.

In Panel 3 you can vary the elements $a$, $b$, $c$ and $d$ of the matrix $\mathbf{A}$ and inspect the behavior of a dynamical system near a fixpoint.

### Memorize this:

It's a good idea to memorize Figure 10 that compiles the different situations:

WLet's apply what we discussed to the following example:

## The Duffing Oscillator

Consider the following dynamical system: [ \dot{x} =y ] [ \dot{y} =-y+x-x^{3} ] This one is known as the Duffing oscillator. It describes the motion of a simple, damped mechanical system that is subject to a nonlinear force. Let's say we have a simple mechanical problem, e.g. a particle of mass $m$ at position $x$ that is subject to a force $F(x)=x-x^{3}$ and a frictional force that is speed dependent $F_{f}=-\gamma v$. We use Newton's second law [ F=ma ] and get [ m\ddot{x}=x-x^{3}-\gamma\dot{x} ] Now this is a one-dimensional system but the dynamical equation involves a second derivative with respect to time. We can also write this as [ \dot{x} =y ] [ m\dot{y} =x-x^{3}-\gamma y ]

identifying the speed $v$ as a second variable. So effectively we have a two-dimensional system of first order. Now if we, for simplicity let $m=1$ and $\gamma=1$ we recover the equations above. We can investigate this system in the usual way, computing the fixpoints, the nullclines plotting the flow etc.. Again, you can try this using the interactive tool. This suggests that we have three fixpoints the trivial one at $(0,0)$ and the two fixpoints at $(\pm1,0).$ The interactive tool suggests that the latter two are stable spirals and the central one is unstable. Let's see if linear stability analysis confirms this observation. We first compute the Jacobian

$$\mathbf{J}(x^{\star},y^{\star})=\left(\begin{array}{cc} 0 & 1\\ 1-3x^{2} & -1 \end{array}\right)$$

Now let's look at the fixpoint $(0,0)$ for which we have

$$\mathbf{J}(0,0)=\left(\begin{array}{cc} 0 & 1\\ 1 & -1 \end{array}\right)$$

which yields $s=-1<0$ and $\Delta=-1<0$. We therefore have an unstable saddle with a positive and negative direction. Now at the other two fixpoints we have

$$\mathbf{J}(\pm1,0)=\left(\begin{array}{cc} 0 & 1\\ -2 & -1 \end{array}\right)$$

which gives $s=-1$ and $\Delta=2>s^{2}/4$ so we have two stable spirals which we expected.

### Asymptotics of the Lotka-Volterra system

Let's do the same for the Lokta-Volterra System:

[ \dot{x} =x(\alpha-\beta y) ] [ \dot{y} =y(\gamma x-\delta) ]

Here we are interested in the stability of the fixpoint $(\delta/\gamma,\alpha/\beta)$. The Jacobian is given by

$$\mathbf{J}(x^{\star},y^{\star})=\left(\begin{array}{cc} \alpha-\beta y & -\beta x\\ \gamma y & \gamma x-\delta \end{array}\right)$$

At the fixpoint we have

$$\mathbf{J}(\delta/\gamma,\alpha/\beta)=\left(\begin{array}{cc} 0 & -\beta\delta/\gamma\\ \gamma\alpha/\beta & 0 \end{array}\right)$$

So the trace is $s=0$ and $\Delta=\alpha\delta$ so we are dealing with a focus. This is interesting. However, unlike some of the other behaviors, this is a behavior that is structurally unstable. Imagine that in addition to the above rhs. of the equations we have a small term

[ \dot{x} =x(\alpha-\beta y)-\epsilon ] [ \dot{y} =y(\gamma x-\delta) ]

This wouldn't change the fixpoint substantially but would imply a negative trace and the focus is replaced by an inward spiral.


text:
# indistinguishable objects

#### schinb65

##### New member
An investor has 20000 to invest among 4 possible investments. Each investment must be a unit of 1000. If the total 20,000 must be invested, how many different investment strategies are possible? What if not all money need to be invested?
I should solve ${\binom{20000+4-1}{20000}}$? I think I need something with the 1000.

Last edited:

#### Jameson

##### Administrator
Staff member
Re: indestinguishable objects

Again, I must preface this with a disclaimer that I'm not confident about my solution.

1) I agree that it's $$\displaystyle \binom{3002}{3000}$$. This isn't that big a number. What calculator are you using?

2) Since these must be in increments of 1000, I think it's really a problem of 20 objects in 4 spaces.

#### schinb65

##### New member
Re: indestinguishable objects

Again, I must preface this with a disclaimer that I'm not confident about my solution.

1) I agree that it's $$\displaystyle \binom{3002}{3000}$$. This isn't that big a number. What calculator are you using?

2) Since these must be in increments of 1000, I think it's really a problem of 20 objects in 4 spaces.
The calculator is a TI 84. I figured out the problem with the calculator. It computes the 3000! first then does the division. This number is too large but if you simplify the numbers it works. Thank you.

#### schinb65

##### New member
Re: indestinguishable objects

Again, I must preface this with a disclaimer that I'm not confident about my solution.

1) I agree that it's $$\displaystyle \binom{3002}{3000}$$. This isn't that big a number. What calculator are you using?

2) Since these must be in increments of 1000, I think it's really a problem of 20 objects in 4 spaces.
Thank you. The 1000 bring the values down to 20 objects gives me the correct answer. Thanks


text:
# Tag Info

18

First, Field strength. This calculation is strictly an electric potential calculation; radiation and induction are safely ignored at 50Hz.* For a 200kV transmission line 20m above ground, the max electric field at ground level is about 1.2 kV/m.** This number is reduced from the naive 200kV/20m=10 kV/m calculation by two effects: 1) The ~1/r variation ...

9

If the power line is 20m high, and has the voltage of 1MV , then the electric field (near ground), very roughly, is on order of 1000/30 kv ~ 30 000 v/m (the numbers are very approximate and the field is complicated because it is a wire near a plate scenario, and wire diameter is unknown but not too small else the air would break down, i.e. spark over, near ...

6

There does not need to be an magnetic field in the inductor for there to be "back emf" (I would prefer "induced emf"). The induced emf is the consequence of a changing magnetic field and not of a magnetic field itself and hence there can be a changing magnetic field even at zero magnetic field (something like a positive acceleration downwards for a ball ...

5

Edit 1: I think I just understood you question: you are actually trying to calculate some sort if “internal” inductance, i.e. the contribution to the inductance of only the field inside the conductor. When calculating the flux, you have to choose a closed path over which you would want the electromotive force, and then integrate the magnetic flux over the ...

4

An induction cooker works with a flat coil setup similar to the one shown below. The diameter of this coil depends on the exact model (in the range of 30 cm). This implies that the field is also extending quite a bit into the air above the cooker. As the large field can be potentially hazardous (think metal ring on finger) there are relatively strict ...

4

Have a look to this old paper: http://www.g3ynh.info/zdocs/refs/NBS/Rosa1908.pdf

4

I see two questions here. The first is why self-inductance is not considered when solving Faraday's law problems, and the second is why an EMF can ever produce a current in a circuit with non-zero self-inductance. I will answer both of these in turn. 1. Why self-inductance is not considered when solving Faraday's law problems Self inductance should be ...

4

In the limit of long times, the currents are steady, so the magnetic fields they create are steady so there is no induced EMF. This situation is usually tagged "steady state". That said, there will be a period of time when you have just switched a circuit on or off during which things have not settled down and then there will in general be effects not seen ...

3

It seems entirely symmetrical. But it isn't symmetrical. For the voltmeter on the right, the 'outer' loop that encloses the magnetic field consists of the voltmeter, leads, and a $100\Omega$ resistor For the voltmeter on the left, the outer loop that encloses the magnetic field consists of the voltmeter, leads, and a $900\Omega$ resistor. Moreover, ...

3

why do we not get a shock? Because the electric resistance of a human body is by orders of magnitude higher than the resistance of the steel pot. why is it that current is converted to heat while it has a good conductor(say, steel) to flow through? According to Maxwell–Faraday equation, changing magnetic field creates the electric field, i.e. the ...

3

Yes you're correct . We willingly don't consider it , it is perfectly correct . However , like in a DC circuit , we think current is established instantaneously . It isn't , first there is a very Induced EMF , and no current flows in circuit . and then there's falling of the EMF with time and it is only at$\infty$ time , that we get the so called current . ...

3

When the switch is opened, the inductor energy is dissipated over time via the internal inductor resistance and the external parallel resistance. The stored energy is radiated away as heat, heating both the inductor and the external parallel resistor. If the external resistance is much larger than the internal inductor resistance, the majority of the ...

3

I believe you simply need 4 op-amps. Here is what you need for differentiation and integration. http://en.wikipedia.org/wiki/Operational_amplifier_applications#Integration_and_differentiation so you simply amplify u by factor 2 with an opamp-circuit then use 2 op-amps for differentiation and integration at the end you should add them all with an opamp ...

3

Given $\theta = 30$ and pure resistance is $20 \Omega$. Using a phase diagram it could be observed that $\tan{\theta} = \frac{\chi_L}{R}$. Find $\chi_L$ and substitute that into $\chi_L = 2 \pi f L$. It's pretty straightforward I believe?

3

Here is an approximation of what you are trying to do. The magnetic field inside a solenoid is given by (from Wikipedia) $B={\mu}nI$ where ${\mu}$ is the permeability of the medium (presumably water in your case), $n$ is the number of turns per unit axial length, and $I$ is the current. The maximum value of $n$ is determined by the thickness of the wire you ...

3

The original question talked about a discrepancy between the result obtained by calculating the flux directly and using the definition $L= \Phi/I$. The confusion arises because of a concept known as "flux linkage". When you calculate the flux enclosed by the region of unit length between r and r+dr, you calculated an expression for flux which you integrated ...

3

Consider an ideal transformer for simplicity (from Wikipedia). (the voltages and currents shown in the picture are phasors) If you consider the above circuit (transformer) as a block in a larger circuit, because it is an ideal transformer you will not have any energy loss in the block. So you can write: $$I_P \times V_P + I_S \times V_S=0$$ and also, ...

3

I know this post is old and has been answered I thought I would post the exact derivation to help out anybody in the future. In order to calculate the the internal inductance of a wire we have to equate the equation for the energy of the magnetic field to the energy from the inductor/inductance. Energy of the $B$ field: $\frac{1}{2\mu}\int B^2dV$, where ...

3

"Mutual inductance" and "self-inductance" are not the same thing, no. Suppose you have two loops (of wire, for example), loop A and loop B. When you run a current through loop A, it produces a magnetic field that fills the space around both loops. The mutual inductance between the two loops, denoted $M_{AB}$, is defined as the ratio of the magnetic flux ...

2

There is a very intuitive way to understand why this must be so. The ideal transformer does not dissipate energy; there is no energy loss (and certainly no energy gain). Thus, if power is delivered (by an external circuit) to the primary, it must be supplied (to another external circuit) by the secondary. If follows that if the primary and secondary ...

2

When calculating the volts per meter of the static field, it's important to assume that the bicycle is conductive (presumably an aluminum frame). Without the bicyclist, one would use image charges to calculate the electric field at the bicycle. The three phases should partially cancel, and Art Brown's calculation seems reasonable, around 1200 volts per ...

2

The relationship $\vec E_c + \vec E_n = 0$ does not hold globally so the two fields are not scalar multiples. Within, and only within, the (ideal) conductor that forms the inductor, the electric field must be zero. Note that the text specifically says so the total electric field ... within the coils must be zero. UPDATE: Consider the general $E$ ...

2

Your magnetic case is similar to the electric case of two opposite charges (balls) of finite density and with no polarization effects (bound charges). At infinite distance you have two equal electrostatic energies $E_1=E_2=E$ (self-energy). At a finite distance $r>2R$ you add to it the energy of electrostatic interaction $U(r)=-Q^2/r$, the latter being ...

2

Books say that ΣΔV=0 around a closed loop KVL holds only if the magnetic flux linking the circuit is unchanging. In ideal circuit theory, it is assumed that circuit elements are ideal lumped elements and the self inductance of the circuit is zero. In other words, we assume that the dimensions of the circuit and circuit elements are arbitrarily small ...

2

Reluctance = $\dfrac{l_e}{\mu A_e}$ where..... mu is the absolute permeability of the material, $\mu_0 \mu_r$ $l_e$ is the circumference of a circle at a radius r and $A_e$ is a small cross sectional area. The circle I refer to only relates to the cross section of the torus and r is the radius from the centre (where the wire is). All these reluctances ...

2

I assume the first part, up to But how exactly does it happen? defines and explains your question, and then you show what you think about it so far? It looks like the point where it goes wrong is about what the inductor does. There is nothing about "split-second" and relativistic, it behaves in a pretty symmetric way to the capacitor. It's "dynamic" ...

1

As the voltage between the capacitor's plates decreases, so should the current flowing through the circuit. I don't follow your reasoning here. Recall that, for an ideal capacitor, we have: $$i_C = C\frac{dv_C}{dt}$$ In words, the current through the capacitor is proportional to the rate of change of the voltage across, not the instantaneous value ...

1

You should write down all the equations of you LC circuit and that may help. The short answer is that when you close the switch and let current flow out of the capacitor, it can't flow right away because the rapidly changing current sets up an opposing voltage in the inductor. V = L di/dt And so the current increases slowly, reaching a maximum (as you ...

1

$\def\ddt{\frac{d}{dt}}\def\l\{\left}\def\r{\right}\def\rmS{S}$The network is over-idealized. Any inductor currents $i_1$, $i_2$ with $i_1+i_2 = \frac{e}{R}$ are possible at steady state. If you are given initial currents $i_1(0),i_2(0)$ then current over time results with $i_\rmS:=i_1+i_2$ from \begin{align} \ddt i_1(t) &= \frac1{L_1}\l(e - R\cdot ...

Only top voted, non community-wiki answers of a minimum length are eligible