As part of an assignment, I have been asked to numerically solve the following 2nd-order differential equation. For those wondering, it is a model of groundwater flow through an aquifer beset on both sides by canals.
$${d^2y(l)\over dl^2} = -{R\over K}y(l)^{-1}$$
$$0 \le l \le L$$
$$y(0) = 3$$
$$y(L) = 2$$
Using Euler's Method. I intend to write a script in MatLab to do this for me. R and K are constants. Since we are talking about a second order (nonlinear?) DE, I started of with rewriting it as a system of first-order differential equations.
Let $x_1(l) = y(l)$. Let $x_2 = x_1'(l) = y'(l)$ Then $x_2'(l) = y"(l)$. The system then becomes:
- $x_1'(l) = x_2(l)$
- $x_2'(l) = -{R\over K}x_1(l)^{-1}$
And this is where I get stuck. All the literature I have read explains the method of rewriting a $n^{th}$-order DE into a system of first-order equations, but only for linear initial value problems. So:
- How do I incorporate the boundary condition at y(L)?
- What to do about the fact I have no boundary conditions for y'(l)? It seems to me that I can never compute the first iteration without knowing a boundary y'(0).
As a quick side-question,
$${d^2y(l)\over dl^2}y(l) = -{R\over K}$$
$${d^2y(l)^2\over dl^2} = -{R\over K}$$
Is just plain wrong, isn't it?
Best Answer
Euler's method is for initial value problems, not boundary value problems. You can try a couple of things.
Shooting method: solve the problem wiyh initial values $y(0)=3$, $y'(0)=a$, where $a$ is a parameter; call $y(x,a)$ the solution. Now you have to solve $y(L,a)=2$. You can do this looking for $a_1,a_2$ such that $y(L,a_1)<3<y(L,a_2)$ and usng the midpoint method.
Finite differences: Take a large $N$, let $h=L/N$ and $l_k=k\,h$, $0\le k\le N$. Let $y_k$ be an approximation of $y(l_k)$. Approximate the differential equation by the system $$ \frac{y_{k+1}-2\,y_k+y_{k-1}}{h^2}=-\frac{R}{K}y_k^{-1},\quad1\le k\le N $$ complemented with $y_0=y(0)=3$, $y_N=y(L)=2$. Solve the system. Since it is a nonlinear system, use Newton-Ralphson.