#Introduction to Numerical Mathematics

On a [0,1] domain find a function u for given source function f and boundary values u_L and u_R such that:

• -u'' = f
• u(0) = u_L
• u(1) = u_R

u'' denotes the second derivative of u

This can be solved purely theoritical but your task is it to solve it numerically on a discretized domain x for N points:

• x = {i/(N-1) | i=0..N-1} or 1-based: {(i-1)/(N-1) | i=1..N}
• h = 1/(N-1) is the spacing

###Input

• f as function or expression or string
• u_L, u_R as floating point values
• N as integer >=2

###Output

• Array, List, some sort of seperated string of u such that u_i == u(x_i)

###Examples

Example 1

Input: f = -2, u_L = u_R = 0, N = 10 (Don't take f=-2 wrong, it is not a value but a constant function that returns -2 for all x. It is like a constant gravity force on our rope.)

There exists an easy exact solution: u = -x*(1-x)

Example 2

Input: f = 10*x, u_L = 0 u_R = 1, N = 15 (Here there is a lot of upwind on the right side)

The exact solution for this states: u = 1/3*(8*x-5*x^3)

Input: f = sin(2*pi*x), u_L = u_R = 1, N = 20 (Someone broke gravity or there is a sort of up- and downwind)

Here the exact solution is u = (sin(2*π*x))/(4*π^2)+1

Input: f = exp(x^2), u_L = u_R = 0, N=30

Note the slight unsymmetry

###FDM

• rewrite -u_i'' = f_i as
• (-u_{i-1} + 2u_i - u{i+1})/h² = f_i which equals
• -u_{i-1} + 2u_i - u{i+1} = h²f_i
• Setup the equations:

• Solve this equation and output the u_i

You may however use any other method to solve the Laplace equation. If you use an iterative method, you should iterate until the residual |b-Au|<1e-6, with b being the right hand side vector u_L,f_1h²,f_2h²,...

###Notes

Depending on your solution method you may not solve the examples exactly to the given solutions. At least for N->infinity the error should approach zero.

###Bonus

###Winning

Introduction to Numerical Mathematics

On a \$[0,1]\$ domain, find a function \$u\$ for given source function \$f\$ and boundary values \$u_L\$ and \$u_R\$ such that:

• \$-u'' = f\$
• \$u(0) = u_L\$
• \$u(1) = u_R\$

\$u''\$ denotes the second derivative of \$u\$

This can be solved purely theoretically but your task is it to solve it numerically on a discretized domain \$x\$ for \$N\$ points:

• \$x = \{\frac i {N-1} : 0 \le i \le N-1\}\$ or 1-based: \$\{\frac {i-1} {N-1} : 0 \le i \le N-1\}\$
• \$h = \frac 1 {N-1}\$ is the spacing

Input

• \$f\$ as a function, expression or string
• \$u_L\$, \$u_R\$ as floating point values
• \$N \ge 2\$ as an integer

Output

• Array, List, some sort of separated string of \$u\$ such that \$u_i = u(x_i)\$

Examples

Example 1

Input: \$f = -2\$, \$u_L = u_R = 0\$, \$N = 10\$ (Don't take \$f=-2\$ wrong, it is not a value but a constant function that returns \$-2\$ for all \$x\$. It is like a constant gravity force on our rope.)

There exists an easy exact solution: \$u = -x(1-x)\$

Example 2

Input: \$f = 10x\$, \$u_L = 0\$, \$u_R = 1\$, \$N = 15\$ (Here there is a lot of upwind on the right side)

The exact solution for this states: \$u = \frac 1 3(8x-5x^3)\$

Input: \$f = \sin(2\pi x)\$, \$u_L = u_R = 1\$, \$N = 20\$ (Someone broke gravity or there is a sort of up- and downwind)

Here the exact solution is \$u = \frac {\sin(2πx)} {4π^2}+1\$

Input: \$f = \exp(x^2)\$, \$u_L = u_R = 0\$, \$N=30\$

Note the slight asymmetry

FDM

• Rewrite \$-u_i'' = f_i\$ as \$\frac {-u_{i-1} + 2u_i - u_{i+1}} {h^2} = f_i\$, which equals \$-u_{i-1} + 2u_i - u_{i+1} = h^2 f_i\$
• Setup the equations:

$$ u_0 = u_L \\ \frac {-u_0 + 2u_1 - u_2} {h^2} = f_1 \\ \frac {-u_1 + 2u_2 - u_3} {h^2} = f_2 \\ \dots = \dots \\ \frac {-u_{n-3} + 2u_{n-2} - u_{n-1}} {h^2} = f_{n-2} \\ u_{n-1} = u_R $$

$$ \begin{pmatrix} 1 & & & & & & \\ -1 & 2 & -1 & & & & \\ & -1 & 2& -1& & & \\ & & \ddots & \ddots & \ddots & & \\ & & & -1 & 2 & -1 & \\ & & & & -1 & 2 & -1 \\ & & & & & & 1 \end{pmatrix} \begin{pmatrix} u_0 \\ u_1 \\ u_2 \\ \vdots \\ u_{n-3} \\ u_{n-2} \\ u_{n-1} \end{pmatrix} = \begin{pmatrix} u_L \\ h^2 f_1 \\ h^2 f_2 \\ \vdots \\ h^2 f_{n-3} \\ h^2 f_{n-2} \\ u_R \end{pmatrix} $$

• Solve this equation and output the \$u_i\$

You may however use any other method to solve the Laplace equation. If you use an iterative method, you should iterate until the residual \$|b-Au| < 10^{-6}\$, with \$b\$ being the right hand side vector \$u_L,f_1 h^2,f_2 h^2, \dots\$

Notes

Depending on your solution method you may not solve the examples exactly to the given solutions. At least for \$N \to \infty\$ the error should approach zero.

Bonus

Winning

You may however use any other method to solve the Laplace equation. If you use an iterative method, you should iterate until the residual |b-Au|<1e-6, with b being the right hand side vector u_L,f_1h²,f_2h²,...

Alternative implementation without Matrix Algebra (using the Jacobi method)

def laplace(f, uL, uR, N):
 h=1./(N-1)
 b=[f(i*h)*h*h for i in range(N)]
 b[0],b[-1]=uL,uR
 u = [0]*N

 def residual():
  return np.sqrt(sum(r*r for r in[b[i] + u[i-1] - 2*u[i] + u[i+1] for i in range(1,N-1)]))

 def jacobi():
  return [uL] + [0.5*(b[i] + u[i-1] + u[i+1]) for i in range(1,N-1)] + [uR]

 while residual() > 1e-6:
  u = jacobi()

 return u

You may however use any other method to solve the Laplace equation. If you use an iterative method, you should iterate until the residual |b-Au|<1e-6, with b being the right hand side vector u_L,f_1h²,f_2h²,...

Note the slight unsymmetry

Note the slight unsymmetry