Commit 73fb02bb authored by Jochen Schulz's avatar Jochen Schulz
Browse files

added exercise 23, missing solutions

parent f66f1a08
%% Cell type:markdown id: tags:
# Exercises 23
## problem 1
Consider the following system of ODE's
$$y_1'(t) = - 0.04 y_1(t) + 10^4 y_2(t)y_3(t)$$
$$y_2'(t) = 0.04 y_1(t) - 10^4 y_2(t) y_3(t) - 3 \cdot 10^7 y_2(t)^2$$
$$y_3'(t) = 3 \cdot 10^7 y_2(t)^2$$
Solve the system with initial value $(1,0,0)$ at $0$ with the time range$0 \leq
t \leq 3$. Use at least two different solvers and compare their results. Try to understand heuristicly why different solvers behave sometimes very differently.
%% Cell type:markdown id: tags:
##problem 2
Implement an explicit Euler method to solve ODE's. The explicit euler method has the parameters $c_i = 0$, $a_{ij} = 0$ and $b_1 = 0$ which then yields the update formula
$$y_{n + 1} = y_n + h f ( t_n , y_ n )$$
Choose a value h for the size of every step and set $t_n = t_0 + n h$ .
%% Cell type:markdown id: tags:
## problem 3
*stiff* ODE's are somewhat *ill-conditioned* in the sense, that many algorithms for solving this type of ODE are then unstable. Consider the initial value problem:
$$y^\prime ( t ) = - 15 y ( t ) , t \ge 0 , y ( 0 ) = 1$$
and solve it with the `odeSolve` and the explicit euler. Look especially at $h= 1/4$ and $h = 1/8$ for the euler method and look what is happening.
%% Cell type:markdown id: tags:
# Exercises 23
## problem 1
Consider the following system of ODE's
$$y_1'(t) = - 0.04 y_1(t) + 10^4 y_2(t)y_3(t)$$
$$y_2'(t) = 0.04 y_1(t) - 10^4 y_2(t) y_3(t) - 3 \cdot 10^7 y_2(t)^2$$
$$y_3'(t) = 3 \cdot 10^7 y_2(t)^2$$
Solve the system with initial value $(1,0,0)$ at $0$ with the time range$0 \leq
t \leq 3$. Use at least two different solvers and compare their results.
t \leq 3$. Use at least two different solvers and compare their results. Try to understand heuristicly why different solvers behave sometimes very differently.
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
##problem 2
Implement an explicit Euler method to solve ODE's. The explicit euler method has the parameters $c_i = 0$, $a_{ij} = 0$ and $b_1 = 0$ which then yields the update formula
$$y_{n + 1} = y_n + h f ( t_n , y_ n )$$
Choose a value h for the size of every step and set $t_n = t_0 + n h$ .
%% Cell type:code id: tags:
``` haskell
import Numeric.Container
import Graphics.Rendering.Chart.Easy hiding (Matrix, Vector, scale)
import Graphics.Rendering.Chart.Backend.Cairo
harmonic w d t v = fromList [y2, -w^2*y1 -2*d*w*y2]
where y1 = v @> 0
y2 = v @> 1
-- For stability: calculate time-step:
--ht = h**2*h**2/( 2*a*(h**2+h**2) )
e_euler :: (Vector Double -> Vector Double) -> Double -> Vector Double -> Vector Double
e_euler f h y = y `add` (h `scale` (f y))
integr :: (Vector Double -> Vector Double) -> Vector Double -> Double -> Int -> [Vector Double]
integr derivative initial_state h n = take n $ iterate (e_euler derivative h) initial_state
```
%% Cell type:code id: tags:
``` haskell
n = 4000
h = 0.005
expEq = integr (harmonic 1 0.1 0) (fromList [1::Double,0]) h n
toRenderable $ do
layout_title .= "Damped harmonic oscillator"
setColors [opaque red]
plot (line "damped (d= 0.1)" [zip (take n [0,0.2..] ) (map (head . toList) expEq)])
```
%% Output
%% Cell type:code id: tags:
``` haskell
-- beware the step size
n = 100
h = 0.2
expEq = integr (harmonic 1 0.1 0) (fromList [1::Double,0]) h n
toRenderable $ do
layout_title .= "Damped harmonic oscillator"
setColors [opaque red]
plot (line "damped (d= 0.1)" [zip (take n [0,0.2..] ) (map (head . toList) expEq)])
```
%% Output
%% Cell type:markdown id: tags:
## problem 3
*stiff* ODE's: Consider the initial value problem:
*stiff* ODE's are somewhat *ill-conditioned* in the sense, that many algorithms for solving this type of ODE are then unstable. Consider the initial value problem:
$$y^\prime ( t ) = - 15 y ( t ) , t \ge 0 , y ( 0 ) = 1$$
and solve it with the odesolve
and solve it with the `odeSolve` and the explicit euler. Look especially at $h= 1/4$ and $h = 1/8$ for the euler method and look what is happening.
%% Cell type:code id: tags:
``` haskell
```
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment