PENG9570 Lecture 5

PENG9570 — Lecture 5 (March 3rd)

Euler’s Method for ODEs

$$y^{\prime }=f(t,y),y\left(t_0\right)=y_0$$

The slope of the solution $y(t)$ is $f(t,y)$.

Linearising on an interval $[t,t+h]$,

$$y(t+h)\approx y(t)+hf(t,y(t)))$$

This is the basis for Euler’s method

Numerical methods: $O\left(h^N\right)$-notation $h$ : Time step Error at time $t=t_n=t_0+nh$ :

$$\left|y_n-y_{\text{Exact }}\left(t_n\right)\right|\le C{h}^N$$

Then the method is said to be $O\left(h^N\right)$ Larger $N$ means higher precision

Numerical Methods for ODEs

$$y^{\prime }=f(t,y),y\left(t_0\right)=y_0$$

Integrating from $t$ to $t_0$, we get the integral form of the IVP: Integral form of IVP:

$$y=y_0+{\int }_{t_0}^{t}f(t,y)dt$$

In general, the integral is evaluated numerically.

Numerical Integration Methods: Riemann Sums

$R_n^{(l)}=\left(f\left(x_0\right)+f\left(x_1\right)+\cdots +f\left(x_{n-1}\right)\right)\Delta x$ $R_n^{(r)}=\left(f\left(x_1\right)+f\left(x_2\right)+\cdots +f\left(x_n\right)\right)\Delta x$

• [?] Left / right Riemann?

Local error $=O\left(h^2\right)$ ⇒ Euler’s method

Euler’s method $y_{n+1}=y_n+hf\left(t_n,y_n\right),y_0=y\left(t_0\right)$ $h$ : Time step Method is $O(h)$ (first order)

Numerical Integration Methods: Midpoint Rule

Mid-point Method

$$M_n=\Delta x\left(f_1^{\ast }+f_2^{\ast }+\cdots +f_{n-1}^{\ast }+f_n^{\ast }\right)$$

Local error $=O\left(h^3\right)$ ⇒ Mid-point method

Estimate using Euler’s method

$$\begin{array}{rl}& {\bar{y}}_{n+1/2}=y_n+\frac{h}{2}f\left(t_n,y_n\right)\\ & y_{n+1}=y_n+hf\left(t_{n+1/2},{\bar{y}}_{n+1/2}\right)\end{array}$$

Method is $O\left(h^2\right)$ (second order)

Numerical Integration Methods: Trapezoidal Rule

$$T_n=\frac{\Delta x}{2}\left(f_0+2f_1+\cdots +2f_{n-1}+f_n\right)$$

Local error $=O\left(h^3\right)$ ⇒ Heun’s method/ Runge-Kutta’s 2. order method

RK2/Heun’s Method

$$\begin{array}{c}{\bar{y}}_{n+1}=y_n+hf\left(t_n,y_n\right),\\ y_{n+1}=y_n+\frac{h}{2}\left(f\left(t_n,y_n\right)+f\left(t_n,{\bar{y}}_{n+1}\right)\right),y_0=y\left(t_0\right)\end{array}$$

Method is $O\left(h^2\right)$ (second order)

Numerical Integration Methods: Simpson’s Rule

$S_{2n}=\frac{\Delta x}{3}\left(f_0+4f_1+2f_2+\cdots +4f_{2n-1}+f_{2n}\right)$ Local error $=O\left(h^5\right)$ ⇒ Runge-Kutta’s 4. order method

Runge-Kutta 4th order Method

$$\begin{array}{c}{k}_1=hf\left(t_n,y_n\right),\\ k_2=hf\left(t_{n+1/2},y_n+\frac{1}{2}{k}_1\right),\\ k_3=hf\left(t_{n+1/2},y_n+\frac{1}{2}{k}_2\right),\\ k_4=hf\left(t_{n+1},y_n+k_3\right),\\ y_{n+1}=y_n+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)\end{array}$$

Method is $O\left(h^4\right)$ (fourth order)

Demonstration of Convergence Properties for the Four Methods

$E$ : Numerical error

$$\begin{array}{c}E\approx C{h}^N\\ \iff \\ \ln E\approx \ln C+N\ln h\end{array}$$

(Linear in $\ln h$ with slope $N$ )

Problem

Consider the IVP

$$y^{\prime }=y^3-\frac{1}{2}y,y(0)=\frac{1}{2}.$$

1. Show that $y(x)=\frac{1}{\sqrt{2e^x+2}}$.
2. Compute numerical solutions on $[0,2]$ for a range of values for $n=$ number of subintervals using (i) Euler’s method, (ii) RK2 (Heun’s method), (iii) the midpoint method and (iv) Runge-Kutta’s fourth order method. |stimate
3. Compute the error $E=\left|y(2)-y_{\text{num }}(2)\right|$ for all methods for all values for $n$ and use these quantities to estimate the convergence rates.

$$y^{\prime }=y^3-\frac{1}{2}y,y(0)=\frac{1}{2}.$$

1. Analytical Solution

Show that $y(x)=\frac{1}{\sqrt{2e^x+2}}$.

Assuming: $y=(2e^x\cdot 2{)}^{-2}$ Then:

$$\begin{array}{rl}& \ast y(0)=\frac{1}{\sqrt{2\cdot e^0+2}}=\frac{1}{\sqrt{4}}=\frac{1}{2}✅\end{array}$$ $$\begin{array}{rl}& \ast y^{\prime }-\left(y^3-\frac{1}{2}y\right)\\ & =-\frac{1}{2}{\left(2e^x+2\right)}^{-3/2}\cdot 2e^x-{\left(2e^x+2\right)}^{-3/2}+\frac{1}{2}{\left(2e^x+2\right)}^{-1/2}\\ & =-e^x{\left(2e^x+2\right)}^{-3/2}-{\left(2e^x+2\right)}^{-3/2}\left(1-\frac{1}{2}\left(2e^x+2\right)\right)\\ & =-e^x{\left(2e^x+2\right)}^{-3/2}-{\left(2e^x+2\right)}^{-3/2}\left(1-e^x-1\right)\\ & =-e^x{\left(2e^x+2\right)}^{-3/2}+e^x{\left(2e^x+2\right)}^{-3/2}=0\end{array}$$

Fill in the transformation from the LHS of the first part of this to its RHS. I'm having some difficulty understanding

ROBOT GENERATED – gemini-3.5-flash

To understand this verification, it helps to first point out a few typographical/OCR errors in your notes that make the math look confusing:

1. First line: $y^{\prime }-\left(y^3-\frac{1}{2}y\right)y$ has an extra $y$ at the end. It should be $y^{\prime }-\left(y^3-\frac{1}{2}y\right)=0$.
2. Second line: $(9e^x+2{)}^{-3/2}$ is a typo for $(2e^x+2{)}^{-3/2}$. Additionally, the exponents in the last two terms should be negative: $-{\left(2e^x+2\right)}^{-3/2}$ and $+\frac{1}{2}{\left(2e^x+2\right)}^{-1/2}$.
3. Fourth line: $-e^x{\left(2e^x+2\right)}^{-7/2}$ is a typo for $-e^x{\left(2e^x+2\right)}^{-3/2}$.

---

Step-by-Step Derivation

Given the analytical solution:

$$y(x)=(2e^x+2{)}^{-1/2}$$

1. Compute the components

Using the chain rule, find the derivative $y^{\prime }$:

$$y^{\prime }=-\frac{1}{2}(2e^x+2{)}^{-3/2}\cdot \frac{d}{dx}(2e^x+2)=-\frac{1}{2}(2e^x+2{)}^{-3/2}\cdot 2e^x=-e^x(2e^x+2{)}^{-3/2}$$

The powers of $y$ are:

$$y^3={\left[(2e^x+2{)}^{-1/2}\right]}^3=(2e^x+2{)}^{-3/2}$$ $$\frac{1}{2}y=\frac{1}{2}(2e^x+2{)}^{-1/2}$$

2. Substitute into the LHS of the ODE ($y^{\prime }-y^3+\frac{1}{2}y$)

$$-e^x(2e^x+2{)}^{-3/2}-(2e^x+2{)}^{-3/2}+\frac{1}{2}(2e^x+2{)}^{-1/2}$$

3. Factor out the common term $(2e^x+2{)}^{-3/2}$

This is the step that usually causes confusion. To simplify the last two terms, factor out $-(2e^x+2{)}^{-3/2}$:

$$-(2e^x+2{)}^{-3/2}+\frac{1}{2}(2e^x+2{)}^{-1/2}=-(2e^x+2{)}^{-3/2}\left[1-\frac{1}{2}(2e^x+2{)}^1\right]$$

(This works because $(2e^x+2{)}^{-3/2}\cdot (2e^x+2{)}^1=(2e^x+2{)}^{-1/2}$.)

4. Simplify inside the brackets

$$1-\frac{1}{2}(2e^x+2)=1-(e^x+1)=-e^x$$

5. Combine and resolve

Now, substitute this simplified expression back into the full equation:

$$-e^x(2e^x+2{)}^{-3/2}-(2e^x+2{)}^{-3/2}(-e^x)$$ $$=-e^x(2e^x+2{)}^{-3/2}+e^x(2e^x+2{)}^{-3/2}=0$$

This confirms that the LHS simplifies to $0$, proving that $y(x)=\frac{1}{\sqrt{2e^x+2}}$ is indeed the correct analytical solution.

Programming: Compute Numerical Solutions

2. Programming: Compute numerical solutions on $[0,2]$ for a range of values for $n=$ number of subintervals using (i) Euler’s method, (ii) RK2 (Heun’s method), (iii) the mid-point method and (iv) Runge-Kutta’s fourth order method.

Exact (curve) Euler (dots) Midpoint (dots) RK2 (dots) RK4 (dots)

3. Programming: Compute the error $E=\left|y(2)-y_{\text{num }}(2)\right|$ for all methods for all values for $n$ and use these quantities to estimate the convergence rates.

Figure 4

Euler, slope=1.0016 RK2, slope=2.0287 Midpoint, slope=2.0276 RK4, slope=4.1342

Example – Model of prey–predator interactions.

Model of prey–predator interactions.

• $x$: prey population,
• $y$: predator population. $x^{\prime }=x(x(1-x)-y)x,y⩾0$ $y^{\prime }=y(x-a),a>0$.

a) Nullclines:

i.e where: $\left(x^{\prime }=0,y^{\prime }=0\right)$

$x^{\prime }=0:\{\begin{array}{l}x(x(1-x)-y)=0\\ \iff x=0\text{ or }y=x(1-x)\end{array}$

$y^{\prime }=0:\{\begin{array}{l}y(x-a)=0\\ \iff y=0\text{ or }x=a.\end{array}$

The case $a>1$: predator death rate is high

(predator death rate is high)

The case $0<a<1$: Equilibrium-ish

Jacobian

How did we get this Jacobian? (oh lol it was already defined, just in reverse order)

ROBOT GENERATED – gemini-3.5-flash

To find the Jacobian matrix $J$, we use the system of differential equations describing the prey-predator interactions given in the note:

$$\begin{array}{rl}{x}^{\prime }& =P(x,y)=x(x(1-x)-y)=x^2-x^3-xy\\ y^{\prime }& =Q(x,y)=y(x-a)=xy-ay\end{array}$$

The Jacobian matrix $J$ is defined as the matrix of first-order partial derivatives of these functions:

$$J=\left[\begin{array}{cc}{P}_x& P_y\\ Q_x& Q_y\end{array}\right]=\left[\begin{array}{cc}\frac{\partial P}{\partial x}& \frac{\partial P}{\partial y}\\ \frac{\partial Q}{\partial x}& \frac{\partial Q}{\partial y}\end{array}\right]$$

Step-by-step Differentiation

1. Top-left entry ($P_x$): Differentiate $P(x,y)$ with respect to $x$ (treating $y$ as a constant):

$$P_x=\frac{\partial }{\partial x}(x^2-x^3-xy)=2x-3x^2-y$$

2. Top-right entry ($P_y$): Differentiate $P(x,y)$ with respect to $y$ (treating $x$ as a constant):

$$P_y=\frac{\partial }{\partial y}(x^2-x^3-xy)=-x$$

3. Bottom-left entry ($Q_x$): Differentiate $Q(x,y)$ with respect to $x$ (treating $y$ as a constant):

$$Q_x=\frac{\partial }{\partial x}(xy-ay)=y$$

4. Bottom-right entry ($Q_y$): Differentiate $Q(x,y)$ with respect to $y$ (treating $x$ as a constant):

$$Q_y=\frac{\partial }{\partial y}(xy-ay)=x-a$$

Final Matrix

Placing these partial derivatives into the Jacobian structure yields the highlighted matrix:

$$J=\left[\begin{array}{cc}2x-3x^2-y& -x\\ y& x-a\end{array}\right]$$

Jacobian:

$$J=\left[\begin{array}{cc}2x-3x^2-y& -x\\ y& x-a\end{array}\right]$$

where, $P(x,y)=x^2(1-x)-xy$ and $Q(x,y)=y(x-a)$,

$$P_x=2x-3x^2-y,P_y=-x,Q_x=y,Q_y=x-a$$

At $(x,y)=(1,0)$:

$J=\left[\begin{array}{cc}2\cdot 1-3\cdot 1^2-0& -1\\ 0& 1-a\end{array}\right]$

$$=\left[\begin{array}{cc}-1& -1\\ 0& 1-a\end{array}\right]$$

Upper triangular: eigenvalues are on the diagonal, ${\lambda }_1=-1,{\lambda }_2=1-a$. $a>1:{\lambda }_1,{\lambda }_2<0$; stable node. $a<1:{\lambda }_1<0,{\lambda }_2>0$; saddle (unstable).

Why are these x,y values selected?

ROBOT GENERATED – gemini-3.5-flash

The values $(x,y)=(1,0)$ are selected because they represent an equilibrium point (or fixed point) of the prey-predator system.

Equilibrium points occur where the system’s rates of change are zero ($x^{\prime }=0$ and $y^{\prime }=0$). We find these by finding the intersections of the system’s nullclines:

1. Prey nullclines ($x^{\prime }=0$): $x=0\text{or}y=x(1-x)$
2. Predator nullclines ($y^{\prime }=0$): $y=0\text{or}x=a$

By intersecting the prey nullcline $y=x(1-x)$ with the predator nullcline $y=0$, we get: $0=x(1-x)⟹x=0\text{or}x=1$

This yields two equilibrium points: $(0,0)$ and $(1,0)$.

Evaluating the Jacobian matrix $J$ at $(1,0)$ allows us to perform linear stability analysis to determine the behavior of the populations near this steady state (e.g., whether it is a stable node or an unstable saddle point depending on the parameter $a$).

At $(x,y)=(a,a(1-a))$ :

$(x,y)=(a,a(1-a))$ (there, $0<a<1$ )

$$\begin{array}{rl}f& =\left[\begin{array}{cc}2a-3a^2-a(1-a)& -a\\ a(1-a)& a-a\end{array}\right]\\ & =\left[\begin{array}{cc}2a-3a^2-a+a^2& -a\\ a(1-a)& 0\end{array}\right]\\ & =\left[\begin{array}{cc}-2a^2+a& -a\\ a(1-a)& 0\end{array}\right]\end{array}$$

Eigenvalues:

$\det (J-\lambda I)=0$

$$\Rightarrow \underset{\text{Characteristic equation }}{\underbrace{{\lambda }^2+\left(2a^2-a\right)\lambda +a^2(1-a)=0}}$$

Define "Characteristic Equation" here

ROBOT GENERATED – gemini-3.5-flash

Based on the provided notes and the highlighted snippet, the characteristic equation is the algebraic equation solved to find the eigenvalues ($\lambda$) of a matrix.

It is defined by setting the determinant of the matrix $(J-\lambda I)$ to zero:

$\det (J-\lambda I)=0$

In the specific context of the prey-predator model evaluated at the equilibrium point $(x,y)=(a,a(1-a))$, expanding this determinant yields the quadratic characteristic equation:

${\lambda }^2+\left(2a^2-a\right)\lambda +a^2(1-a)=0$

where:

• $J$ is the Jacobian matrix.
• $I$ is the identity matrix.
• $\lambda$ represents the eigenvalues.

$$\Rightarrow \lambda =\frac{a}{2}\left[-(2a-1)\pm \sqrt{(2a-1{)}^2-4(1-a)}\right]$$

Result:

• Stable node: For $\frac{\sqrt{3}}{2}<a<1$, we have two real eigenvalues, both negative.
• Stable spiral (sink): For $\frac{1}{2}<a<\frac{\sqrt{3}}{2}$, we have complex conjugate eigenvalues with negative real part.
• Unstable spiral (source): For $0<a<\frac{1}{2}$, we have complex conjugate eigenvalues with positive real part.

Where is the bound on the a values coming from? Is it from solving for a with the eigenvalues?

ROBOT GENERATED – gemini-3.5-flash

Yes, the bounds on the $a$ values come directly from solving for and analyzing the eigenvalues ($\lambda$).

From the characteristic equation in the note: $\lambda =\frac{a}{2}\left[-(2a-1)\pm \sqrt{(2a-1{)}^2-4(1-a)}\right]$

By simplifying the discriminant (the term inside the square root), we get: $\Delta =(2a-1{)}^2-4(1-a)=4a^2-4a+1-4+4a=4a^2-3$

1. Real vs. Complex Eigenvalues (The $\frac{\sqrt{3}}{2}$ bound)

• Real eigenvalues occur when the discriminant is positive ($\Delta >0$): $4a^2-3>0⟹a^2>\frac{3}{4}⟹a>\frac{\sqrt{3}}{2}$
• Complex conjugate eigenvalues occur when the discriminant is negative ($\Delta <0$): $4a^2-3<0⟹a^2<\frac{3}{4}⟹a<\frac{\sqrt{3}}{2}$

2. Stable vs. Unstable Spirals (The $\frac{1}{2}$ bound)

When the eigenvalues are complex ($a<\frac{\sqrt{3}}{2}$), their real part is determined by the term outside the square root: $\text{Re}(\lambda )=-\frac{a}{2}(2a-1)$

Since $a>0$, the sign of the real part depends entirely on the term $-(2a-1)$:

• Negative real part (Stable spiral / Sink): $-(2a-1)<0⟹2a-1>0⟹a>\frac{1}{2}$
• Positive real part (Unstable spiral / Source): $-(2a-1)>0⟹2a-1<0⟹a<\frac{1}{2}$

Combining these conditions yields the exact bounds presented in the results.

Result Illustrations

$$\frac{\sqrt{3}}{2}<a<1$$

$$\frac{1}{2}<a<\frac{\sqrt{3}}{2}$$

$$0<a<\frac{1}{2}$$

Hopf Bifurcation Diagram

So if given:

$$\begin{array}{rl}& x^{\prime }=x(x(1-x)-y)\\ & y^{\prime }=y(x-a)\end{array}$$

Then, Hopf bifurcation at $a=0.45$ (he writes 0.5 in slide). Equilibrium becomes unstable and a limit cycle emerges.

Phase Portrait

Blue: $a=1.2$ Yellow: $a=0.85$ Green: $a=0.6$ Red: $a=0.4$

Time-Series Plots

Phase plots using different methods

Euler’s method predicts wrong closed curve solution.

Partial Differential Equations (PDES)

$Functionsofonevariable:u(x)$ $Functionsofmorethanonevariable:\{\begin{array}{l}u(x,y)\\ u(x,t)\\ u(x,y,t)\\ u(x,y,z,t)\end{array}$

Partial Derivatives

$u^{\prime }(x)=\frac{du}{dx}$ : slope of tangent of $u$ at $x$ $u=u(x,y):$

• $u_x=\frac{\partial u}{\partial x}$ : slope of tangent in the $x$-direction at $(x,y)$ ($y$ fixed)
• $u_y=\frac{\partial u}{\partial y}$ : slope of tangent in the $y$-direction at $(x,y)$ ($x$ fixed)

Notation for second order derivatives:

$$\begin{array}{rl}& u_{xx}=\frac{{\partial }^{2}u}{\partial x^2}\\ & u_{yy}=\frac{{\partial }^{2}u}{\partial y^2}\\ & u_{xy}=\frac{{\partial }^{2}u}{\partial x\partial y}\\ & u_{yx}=\frac{{\partial }^{2}u}{\partial y\partial x}\end{array}$$

Example 1 – Heat equation

Heat equation. $K$ : heat conductivity

$$u_t=k{u}_{xx}u=u(x,t)$$

Example 2 – Fisher’s equation

Fisher’s equation. $D$: diffusion constant

$$u_t=D{u}_{xx}+ru\left(1-\frac{u}{k}\right)$$ $$u=u(x,t)$$

Example 3 – Solve: $u_{tx}=tx$

Solve: $u_{tx}=tx,u=u(x,t)$

$$\left[u_{tx}={\left(u_t\right)}_x={\left(u_x\right)}_t\right]$$

We write the PDE as

$${\left(u_t\right)}_x=tx$$

We introduce $v=u_t$ :

$$\underset{\begin{array}{c}\text{ Can be considered to be an ODE with parameter t}\end{array}}{\underbrace{v_x=tx}}$$ $$\begin{array}{rl}\Rightarrow v& =\int txdx=t\int xdx\\ & =t\cdot \frac{1}{2}{x}^2+\underset{\begin{array}{c}\text{ Arbitrary }\\ \text{ function of t(?). }\end{array}}{\underbrace{C(t)}}\end{array}$$

We have

$$\underset{\text{ODE with parameter }x.}{\underbrace{u_t=\frac{1}{2}t{x}^2+C(t)}}$$ $$\begin{array}{rl}\Rightarrow u& =\int \left(\frac{1}{2}t{x}^2+c(t)\right)dt\\ & =\frac{1}{2}{x}^2\int tdt+\int c(t)dt\\ & =\frac{1}{2}{x}^2\cdot \frac{1}{2}{t}^2+\underset{\bar{c}(t)}{\underbrace{\int c(t)dt}}+D(x)\\ & =\underset{GeneralSolution}{\underbrace{\frac{1}{4}{x}^2{t}^2}}+\underset{Arbitrary}{\underbrace{\bar{c}(t)+D(x)}}\end{array}$$

What does the bar over c(t) signify? Why do we add arbitrary functions?

ROBOT GENERATED – gemini-3.5-flash

Based on the provided notes and highlighted mathematical steps, here are the answers to your questions:

1. What does the bar over $c(t)$ signify?

The bar over $c(t)$, written as $\bar{c}(t)$, signifies the antiderivative (or integral) of the arbitrary function $c(t)$ with respect to $t$.

As shown in the snippet: $\bar{c}(t)=\int c(t)dt$

Because $c(t)$ is an arbitrary function of $t$, its integral $\bar{c}(t)$ is also just another arbitrary function of $t$. The bar notation is used to simplify and clean up the final expression of the general solution.

---

2. Why do we add arbitrary functions?

In partial differential equations (PDEs), when you integrate with respect to one variable, you must treat all other independent variables as constants:

1. Integrating with respect to $x$: Since $t$ is treated as a constant, the “constant of integration” does not have to be a single number; it can be any function that depends only on $t$, denoted as $C(t)$. This is because its partial derivative with respect to $x$ is zero: $\frac{\partial }{\partial x}[C(t)]=0$

2. Integrating with respect to $t$: Since $x$ is treated as a constant, the constant of integration can be any arbitrary function of $x$, denoted as $D(x)$, because: $\frac{\partial }{\partial t}[D(x)]=0$

Adding these arbitrary functions is necessary to capture the complete, general solution of the PDE.

Example 4 – Solve : $u_{tx}+2u_x=0$

$$\begin{array}{rl}{u}_{tx}+2u_x& =0,u=u(x,t)\\ {\left(u_t\right)}_x+2u_x& =0\\ {\left(u_t+2u\right)}_x& =0\end{array}$$ $$\begin{array}{r}\Rightarrow \underset{\text{Possible to solve but tedious}}{\underbrace{u_t+2u=\underset{\begin{array}{c}\text{ Arbitrary }\\ \text{ function }\end{array}}{\underbrace{F(t)}}}}\end{array}$$

Possible to solve, but tedious.

Introduce $v=u_x$ :

$$\begin{array}{rl}& u_{tx}+2u_x=0\\ & {\left(u_x\right)}_t+2u_x=0\\ & \underset{\text{ODE with parameter x}}{\underbrace{v_t+2v=0}}\end{array}$$

Then:

$$\begin{array}{rl}& v_t=-2v\\ & v=C(x)e^{-2t}\\ & \underset{\text{ODE with parameter t}}{\underbrace{u_x=C(x)e^{-2t}}}\end{array}$$

• [?] How does it go from step 1 to 2? We derive wrt to $t$ but then where does the whole RHS come from?

Continuing:

$$\begin{array}{rl}\Rightarrow u& =\int C(x)e^{-2t}dx\\ & =e^{-2t}\underset{\bar{C}(x)}{\underbrace{\int C(x)dx}}+D(t)\\ & =\bar{C}(x)e^{-2t}+D(t)\end{array}$$

General solution: $\bar{C}(x)e^{-2t}+D(t)$

Arbitrary functions: $\bar{C}(x),D(t)$

Conservation laws

Note

I’m not exactly sure what is being asked for here

$\begin{array}{l}U(x,t)\text{ : density of substance }\\ \text{ A: cross-sectional area }\end{array}$

I. Amount of substance in $[a,b]$ at time $t$

$$={\int }_a^{b}u(x,t)\cdot Adx$$

II. Net rate of substance flowing into $[a,b]$

$$\begin{array}{rl}& =A\cdot f(a,t)-A\cdot f(b,t)\\ & =A(f(a,t)-f(b,t))\\ & =-A{\int }_a^b{f}_x(x,t)dx\end{array}$$

III. Rate of production of substance in $[a,b]$

$$={\int }_a^b\underset{\begin{array}{c}\text{ Production }\\ \text{ rate at }\end{array}(x,t)}{\underbrace{f(x,t,u)}}Adx$$

IV. Change of amount of substance per time unit

$$\begin{array}{rl}& =\frac{d}{dt}{\int }_a^{b}u(x,t)Adx\\ & ={\int }_a^b{u}_t(x,t)Adx\end{array}$$

Conservation Law:

Net rate of substance into $[a,b]$ + rate of production inside $[a,b]$ = change of amount in $[a,b]$

$$i.e.:\mathbb{II}+\mathbb{II}=\mathbb{IV}$$

Mathematically:

$$\begin{array}{rl}& -A\cdot {\int }_a^b{f}_x(x,t)dx\\ & +{\int }_a^{b}f(x,t,u)Adx\\ & ={\int }_a^b{u}_t(x,t)Adx\end{array}$$ $$\begin{array}{rl}& \Rightarrow A\left[-{\int }_a^b{y}_{x}dx+{\int }_a^{b}fdx-{\int }_a^b{u}_{t}dx\right]=0\\ & \Rightarrow {\int }_a^b\left(u_t(x,t)+y_x(x,t)-f(x,t,u)\right)dx=0\\ & \text{Conservation law on integral form }\end{array}$$

The interval $[a,b]$ is chosen arbitrarily ⇒ the integrand must be Zero!

$$\begin{array}{rl}& u_t+f_x-f=0\\ & u_t+f_x=f\\ & \text{I.e. the Conservation law on differential form.}\end{array}$$

Fick’s law:

Note

I think I have a topics Fick’s law

• Fix this section!

Fick’s law is as follows:

$$J=-D{u}_x$$

$$\begin{array}{rl}{J}_x& ={\left(-D{u}_x\right)}_x\\ & =-D{u}_{xx}\end{array}$$

Insertion:

$$\begin{array}{rl}& u_t-D{u}_{xx}=f\\ & \Rightarrow u_t=D{u}_{xx}+f\\ & \text{ Diffusion PDE (with source term) }\end{array}$$

• [?] I’m not certain what is happening here

Boundary and Initial Conditions:

The initial boundary Value problem (IBVP) Consider:

$$u_t=k{u}_{xx},0<x<l,t>0.$$

To determine the solution Uniquely we need:

• an initial function

$$u(x,0)=f(x),0<x<l.$$

• boundary conditions

$$\begin{array}{rl}& u(0,t)=g(t),t>0\\ & u(l,t)=h(t)\end{array}$$

Dirichlet Conditions:

$f,g$ and $h$ are prescribed functions (as we have here).

Insulated boundary conditions: $u_x(0,t)=u_x(l,t)=0\text{(zero-flux)}$ Neumann boundary conditions:

$$\begin{array}{r}{u}_x(0,t)=g(t),\\ u_x(l,t)=h(t)\end{array}$$

An explicit scheme for the heat / diffusion PDE / IBVP

IBVP: $u_t=u_{xx},t>0,0<x<1$ IC: $u(x,0)=g(x),0<x<1$. BC: $\{\begin{array}{l}u(0,t)=0\\ u(1,t)=0\end{array}$ $[0,1]$ is separated into $n+1$ subintervals of width $\Delta x=\frac{1}{n+1}$. $n=$ number of interior points for one value of $t$. $x_i=i\Delta x$ $t_j=j\Delta t$

Approximations of $u_t$ and $u_{xx}$

$$u_t(x,t)=\underset{\Delta t\to 0}{\lim }\frac{u(x,t+\Delta t)-u(x,t)}{\Delta t}$$

If $\Delta t$ is small:

$$\begin{array}{rl}{u}_t& \approx \frac{u(x,t+\Delta t)-u\left(x_1,t\right)}{\Delta t}\\ t_{t=tj}^{=ti}-t_{j+1}& =\frac{u\left(x_i,t_j+\Delta t\right)-u\left(x_i,t_j\right)}{\Delta t}\\ & =\frac{u\left(x_i,t_{j+1}\right)-u\left(x_i,t_j\right)}{\Delta t}\end{array}$$

We introduce $V_{i,j}$ as an approximation to $u\left(x_i,t_j\right)$ j

$$v_{i,j}\approx u\left(x_i,t_j\right)$$ $$\Rightarrow u_t\left(x_i,t_j\right)\approx \frac{v_{i,j+1}-v_{i,j}}{\Delta t}$$

Three-point approximation for ${\omega }^{\prime \prime }(x)$ :

$$w^{\prime \prime }(x)\approx \frac{w(x+\Delta x)-2w(x)+w(x-\Delta x)}{\Delta x^2}$$

Thus, $U_{xx}\left(x_i,t_j\right)x_{i+1}$ $\approx \frac{u\left({\tilde{x}}_i+\Delta x_j,t_j\right)-2u\left(x_i,t_j\right)+u\left(\frac{\prime \prime }{x_i-\Delta x,t_j}\right)}{\Delta x^2}$ $=\frac{u\left(x_{i+1},t_j\right)-2u\left(x_i,t_j\right)+u\left(x_{i-1},t_j\right)}{\Delta x^2}$

$$\approx \frac{V_{i+1,j}-2V_{i,j}+V_{i-1,j}}{\Delta x^2}$$

Since $u_t=u_{xx}$ :

$$\begin{array}{rl}& \frac{v_{i,j+1}-v_{i,j}}{\Delta t}\\ & =\frac{v_{i+1,j}-2v_{i,j}+v_{i-1,j}}{\Delta x^2}\end{array}$$

Discretized Version of

$$u_x=u_{xx}$$

Derivation of a matrixvector formulation

Let $\alpha =\frac{\Delta t}{\Delta x^2}$

$$\frac{V_{i,j+1}-V_{i,j}}{\Delta t}=\frac{V_{i+1,j}-2V_{i,j}+V_{i-1,j}}{\Delta x^2}$$

The explicit scheme For $i=1$ (application of $B(1;x=0$ )

$$\begin{array}{rl}& V_{1,j+1}=V_{1,j}+\alpha \left(V_{2,j}-2v_{1,j}+v_{0,j}\right)\\ & V_{0,j}\approx u\left(x_0,t_j\right)=u\left(0,t_j\right)=0\end{array}$$ $$\begin{array}{rl}\Rightarrow v_{1,j+1}& =v_{1,j}+\alpha \left(v_{2,j}-2v_{1,j}\right)\\ & =(1-2\alpha )v_{1,j}+\alpha v_{2,j}\end{array}$$

For $i=n$ (application of $B(2:x=1)$

$$\begin{array}{rl}{v}_{n,j+1}& =v_{n,j}+\alpha \left(v_{n+1,j}-2v_{n,j}+v_{n-1,j}\right)\\ v_{n+1,j}& \approx u\left(x_{n+1},+j\right)=u(1,t,)=0)\\ v_{n,j+1}& =v_{n,j}+\alpha \left(-2v_{n,j}+v_{n-1,j}\right)\\ & =\alpha v_{n-1,j}+(1-2\alpha )v_{n,j}\end{array}$$

General case:

$$V_{i,j+1}=\alpha V_{i-1,j}+(1-2\alpha )V_{i,j}+\alpha V_{i+1,j}$$

The equations read, structurally:

$$\begin{array}{rl}& V_{1ij+1}=(1-2\alpha )V_{1ij}+\alpha V_{2ij}\\ & V_{2ij+1}=\alpha V_{1ij}+(1-2\alpha )V_{ij}+\alpha V_{3ij}\end{array}$$ $$\begin{array}{lr}{V}_{n-1;1+1}=& \alpha V_{n-2,j}+(1-2\alpha )V_{n-1,j}+\alpha V_{V_{ij}}\\ V_{n-1j}=& \alpha V_{n-1,j}+(1-2\alpha {)}_{nij}\end{array}$$

Linear system of equations:

$$\left[\begin{array}{c}{v}_{1,j+1}\\ v_{2,j+1}\\ ⋮\\ v_{u_{1,j+1}}\end{array}\right]=\left[\begin{array}{ccccc}1-2\alpha & \alpha & & & \\ \alpha & 1-2\alpha & \alpha & & \\ & \alpha & \ddots & \ddots & \\ & & \ddots & \ddots & \alpha \\ & & & \alpha & 1-2\alpha \end{array}\right]\left(\begin{array}{c}{v}_{1,j}\\ v_{2,j}\\ ⋮\\ v_{u_{1,j}}\end{array}\right)$$

Define ${\vec{V}}_j=\left[\begin{array}{c}{V}_{1,j}\\ V_{2,j}\\ ⋮\\ V_{n,j}\end{array}\right]$ $\hat{A}=\left[\begin{array}{cccc}1-2\alpha & \alpha & & 0\\ \alpha & 1-2\alpha & \ddots & \\ & \ddots & \ddots & \alpha \\ 0& & & \alpha \end{array}\right]$ The explicit scheme:

$${\vec{V}}_{j+1}=\hat{A}{\vec{V}}_j$$ $$\begin{array}{rl}& \text{ So, }\\ & {\vec{V}}_0=\underset{k_{\text{nown }}}{\underbrace{\left[\begin{array}{c}{V}_{1,0}\\ V_{2,0}\\ ⋮\\ V_{n,0}\end{array}\right]=\left[\begin{array}{c}g\left(x_1\right)\\ g\left(x_2\right)\\ ⋮\\ g\left(x_n\right)\end{array}\right]}}\\ & {\vec{V}}_1=\hat{A}{\vec{V}}_0\\ & {\vec{V}}_2=\hat{A}{\vec{V}}_1=\hat{A}\hat{A}{V}_0={\hat{A}}^2{V}_0\\ & {\vec{V}}_j={\hat{A}}^j{\vec{V}}_0\text{("Solution formula")}\end{array}$$