﻿ Euler Math Toolbox - Examples

# Singular Boundary Value Problems

This notebook is based on a worksheet by Radovan Omorjan.

We will solve the boundary value problem for the second order ordinary differential equation given in the form

```  y" + g1(x,y)*y' + g2(x,y)*y = g3(x)
```

with the boundary conditions

```  A*y(x=a) + B*y'(x=a) = C
```
```  A1*y(x=b) + B1*y'(x=b) = C1
```

# Example

```  y" - phi*y^n = 0
```

Boundary conditions

```  y(x=0)= 0.1
y(x=1)= 1
```

phi,n are given constants

Defining the function of the derivatives, we set y=y, y=y' therefore

```  y'=y
y'=g3(x)-g2(x,y)*y-g1(x,y)*y
```
```>function deriv(x,y) := [y,g3(x)-g2(x,y)*y-g1(x,y)*y]
```

The user has to define functions g1(x,y),g2(x,y),g3(x), f0(x) and the guess value for y'(0).

```>function g1(x,y) := 0
>function g2(x,y) := -phi*y^n/y
>function g3(x) := 0
```

Defining the function in single variable which has zero values in the solution (see the given boundary conditions)

```>function f0(ydguess) := lsoda("deriv",[0,1],[0.1,ydguess])[1,2]
```

We use this values.

```>phi=4; n=1;
```

Let us try to values.

```>f0(0), f0(1)
```
```0.376219569124
2.18964977308
```
```>res=secant("f0",0,2,y=1)
```
```0.343978185375
```

Should be close to 1.

```>f0(res)
```
```1
```

Plotting the result!

```>t=0:0.1:1; y=lsoda("deriv",t,[0.1,res]); ...
plot2d(t,y,color=10,points=1,style="o"); ...
insimg(25,"");
``` This simple example of boundary value problem for second order ODE could also be solved in Maxima analytically.

```>function g1(x,y) &= 0
```
```                                  0

```
```>function g2(x,y) &= -phi*y^n/y
```
```                                    n - 1
- phi y

```
```>function g3(x) &= 0
```
```                                  0

```
```>equ &= 'diff(y,x,2) + g1(x,y)*'diff(y,x) + g2(x,y)*y = g3(x)
```
```                            2
d y        n
--- - phi y  = 0
2
dx

```
```>res &= ode2(subst([n=1,phi=4],equ),y,x)
```
```                               2 x        - 2 x
y = %k1 E    + %k2 E

```
```>res1 &= rhs(bc2(res,x=0,y=1/10,x=1,y=1))
```
```                     2       2 x     4       2   - 2 x
(10 E  - 1) E      (E  - 10 E ) E
---------------- + -------------------
4                   4
10 E  - 10          10 E  - 10

```

Plotting numerical (points) and analytical solutions (line)

```>x=0:0.01:1; plot2d(x,res1(x),add=true,color=13); ...
insimg(25,"");
``` # Another problem

We solve a very special case, where the derivative function has a singularity at the boundary.

```  y'' = 2y'/x - 4y
y(0)=0
y(1)=1
```
```>function f(x,y) := [ y, 2*y/x-4*y ]
```

The singularity is removed if y'(0)=0. So this function might have a nice solution. However, y(0)=0, and y'(0)=0 always yield the zero function for normal problems (e.g. Lipschitz continous f). Thus we expect problems.

In fact there is no proper way to shoot from 0 to 1. So we do a backward shooting.

Setup the shooting from x=1, y=1, back to 0, a function depending on y'(1).

```>function g(der) := lsoda("f",[1,0],[1,der])[1,-1]
>g(1)
```
```0.236949825926
```

Find a derivative y'(1)=der with y(0)=0.

```>der1 = secant("g",2.5)
```
```2.0884291999
```

Plot the correct solution.

```>x=1:-0.01:0; y=lsoda("f",x,[1,der1]); plot2d(x,y):
``` Plot two shootings. All of them are defined across x=0.

```>x=1:-0.01:-1; y=lsoda("f",x,[1,3]); plot2d(x,y); ...
insimg;
``` # Example

Let us try a simpler example.

```  y'' = 2*y'/x-4
y(0) = 0
y(1) = 1
```
```>function f(x,y) := [ y, 2*y/x-4 ]
>function g(der) := lsoda("f",[1,0],[1,der])[1,-1]
>der1 = secant("g",2.5)
```
```1.00000000001
```

Plot the correct solution.

```>x=1:-0.01:0; y=lsoda("f",x,[1,der1]); plot2d(x,y):
``` We can solve this with Maxima.

```>sol &= ode2('diff(y,x,2)=2*'diff(y,x)/x-4,y,x)
```
```                                3      2   %k1
y = %k2 x  + 2 x  - ---
3

```
```>function y(x,%k1,%k2) &= rhs(sol)
```
```                              3      2   %k1
%k2 x  + 2 x  - ---
3

```
```>&diff(y(x,a,b),x,2) = expand(-2*diff(y(x,a,b),x)/x-4)
```
```                       6 b x + 4 = - 6 b x - 12

```
```>&solve(y(1,a,b)=1,[a,b]);  ...
function y1(x,c) &= y(x,a,b) with % with %rnum_list=c
```
```                          3      2   - 3 c - 3
c x  + 2 x  + ---------
3

```
```>&expand(y1(1,c))
```
```                                  1

```
```>&diff(y1(x,c),x) with x=1
```
```                               3 c + 4

```
```>plot2d("y1(x,-4/3)",0,1):
``` # Example

The following problem does not have a solution.

```  y'' = -2*y'/x-4
y(0) = 0
y(1) = 1
```

From the general solution you see that the boundary conditions cannot hold.

```>sol &= ode2('diff(y,x,2)=-2*'diff(y,x)/x-4,y,x)
```
```                                 2
2 x    %k2
y = - ---- + --- + %k1
3      x

```

# Another Example

The following example has the solution y(x)=sqrt(x)

```  y'' = -1/2*y'/x
y(0)=0
y(1)=1
```
```>&sqrt(x); &diff(%,x,2) + 1/2 *diff(%,x)/x
```
```                                  0

```

The shooting method is hard to do, since lsoda will fail, if x is too close to 0.

```>function f(x,y) := [ y, -1/2*y/x ]
```

Let us try two shootings from x=1, but only to 0.01

```>x=1:-0.01:0.01; ...
y1=lsoda("f",x,[1,1/4]); ...
y2=lsoda("f",x,[1,3/4]); ...
plot2d(x,y1_y2):
``` ```>function h(der) := lsoda("f",[1,1e-4],[1,der])[1,2];
```

That is almost the best we can do. The correct answer is 0.5.

```>der0:=bisect("h",1/4,3/4)
```
```0.505050504939
```

# Yet another problem

The problem changes again, if we add another term to the solution.

```>function f(x,y) := [ y, -2*y/x + 4*y ]
```

Let us try some shooting.

```>x=1:-0.01:0.01; ...
y1=lsoda("f",x,[1,1.2]); ...
y2=lsoda("f",x,[1,1]); ...
plot2d(x,y1_y2):
``` ```>function h(der) := lsoda("f",[1,0.01],[1,der])[1,2];
>der0:=bisect("h",1,1.2)
```
```1.07773431527
```

The solution with y(0.01)=0 looks like this.

```>y=lsoda("f",x,[1,der0]); ...
plot2d(x,y):
``` There is a nice solution to this equation.

```>function f(x,a,b) &= a*exp(-2*x)/x + b*exp(2*x)/x
```
```                             2 x      - 2 x
b E      a E
------ + --------
x         x

```

We can solve it for y(1)=1, y(s)=0.

```>sol &= solve([f(1,a,b)=1,f(s,a,b)=0],[a,b])
```
```                           4 s + 2           2
E                 E
[[a = - ---------, b = ---------]]
4    4 s       4    4 s
E  - E         E  - E

```

And make the solution a function.

```>function fsol(x,s) &= factor(f(x,a,b) with sol)
```
```                2 - 2 x   x    s    x    s    2 x    2 s
E        (E  - E ) (E  + E ) (E    + E   )
- ------------------------------------------
s        s        2 s    2
(E  - E) (E  + E) (E    + E ) x

```

Let us check.

```>&ratsimp(diff(fsol(x,s),x,2)+2*diff(fsol(x,s),x)/x-4*fsol(x,s))
```
```                                  0

```

Add the solution to the plot.

```>plot2d("fsol(x,0.01)",color=red,add=true):
``` The limit for s->0 looks quite different.

```>plot2d("fsol(x,0)",0,1):
``` Here, we have the function.

```>&ratsimp(fsol(x,0))
```
```                         - 2 x   4 x + 2    2
E      (E        - E )
----------------------
4
(E  - 1) x

```

We can also solve it for y'(s)=0. Let us define the function for the derivative.

```>function fd(x,a,b) &= diff(f(x,a,b),x)
```
```                   2 x      2 x        - 2 x      - 2 x
2 b E      b E      2 a E        a E
-------- - ------ - ---------- - --------
x          2         x            2
x                      x

```

And solve it.

```>sol &= solve([f(1,a,b)=1,fd(s,a,b)=0],[a,b])
```
```                              4 s + 2
(2 s - 1) E
[[a = -----------------------------,
4 s    4
(2 s - 1) E    + E  (2 s + 1)
2
E  (2 s + 1)
b = -----------------------------]]
4 s    4
(2 s - 1) E    + E  (2 s + 1)

```

Make the solution a function.

```>function fsol(x,s) &= factor(f(x,a,b) with sol)
```
```              2 - 2 x       4 x    4 x        4 s    4 s
E        (2 s E    + E    + 2 s E    - E   )
--------------------------------------------
4 s    4 s      4      4
(2 s E    - E    + 2 E  s + E ) x

```
```>plot2d("fsol(x,0.1)",0,1):
``` The solution for y'(0)=0 happens to be the same as y(0)=0.

```>&fsol(x,0)
```
```                          2 - 2 x   4 x
E        (E    - 1)
-------------------
4
(E  - 1) x

```

Examples