<content> DSolRes </content>

affa75a0-f5fe-4c3b-b08b-9592108a8f2b 401 4 5 false false 0 decimal

&[DATE] &[TIME] - &[FILENAME]

Solving ODE's with the theta or the trapezoidal rule

The theta method uses a parameter 0 < θ < 1 for solve ODE's. Special values of θ are:

θ = 0: Backward Euler
θ = 1/2: Trapezoid method
θ = 1: Forward Euler

<content> DSolRes </content>

θ rule + Newton Raphson ODE Solver Ds toe xi θ N DSolRes m xi length : r 1 m range : c r : hId h m identity * eval : 14 mat X xi r el : T toe min 11 mat eval : dt toe max T 1 el - N / eval : J 0 : 14 mat k 1 N range to T k el : xo X k col : 12 mat t to dt + : x xo dt to xo D * + : 12 mat x res x xo - t to - (θ Ds to xo feval * 1 θ - (Ds t x feval * + (* - : R x res : iter 1 MaxI range R normi εy ≤ break J r c el 1 h / x hId c col + res r el R r el - (* eval : Δx J 1 -^R * - eval : Δx normi εx ≤ break x x Δx + eval : R x res eval : 12 mat if 31 line if for T k 1 + el t : X c k 1 + el x c el : 12 mat 61 line for T X transpose augment 41 line :

The same, but using alg lib solver: faster.

θ rule + Newton Raphson ODE Solver Ds toe xi θ N DSolRes m xi length : r 1 m range : c r : hId h m identity * eval : 14 mat X xi r el : T toe min 11 mat eval : dt toe max T 1 el - N / eval : 13 mat k 1 N range to T k el : xo X k col : 12 mat t to dt + : x xo dt to xo D * + : 12 mat x res x xo - t to - (θ Ds to xo feval * 1 θ - (Ds t x feval * + (* - : T k 1 + el t : X c k 1 + el x res al_nleqsolve c el : 12 mat 41 line for T X transpose augment 41 line :

<content> Utils </content>

f# x# feval f# num2str (x# num2str) concat str2num 11 line :

f# x# y# feval f# num2str (x# num2str, y# num2str) concat str2num 11 line :

M MCols C 0 : k 1 M cols range C 1 k el M k col : for C 31 line :

Defaults

h 1012 -^:

εy 1014 -^:

εx 1014 -^:

0107 -^:

MaxI 100 :

<content> DSolRes Example </content>

Example: Exact linear ODE

x'' A' x * A x' * + B' + ≡

Given

to 0 :

te 3 :

xo 0112 mat :

N 15 :

t B e 0.2 t * -^t * :

t A 0.1 t * :

let

t A' t A t diff : 110 /

t IA τ2 A τ2 to t int :

t B' t B t diff : 5 t - 5 e 5 nthroot t^* /

Convert the ODE into
a system

t x D x 2 el t A' x 1 el * t A x 2 el * + t B' + 21 mat :

Exact solution, numerically evaluated

x.s 10 B - τ B + (e τ IA -^* τ to t int e t IA^* :

t x.s x.s :

The numerical evaluation for the integrals could be very slow. Enable this for speed up the calculations, only if Maple can find a close solution and SMath recognize all functions.

t x.s maple x.s feval : 5 sqrt 25 sqrt 1 t 4 t + (* 20 / - exp - (* π sqrt 5 sqrt 2 t + (* 10 / erf - 15 sqrt / erf + (* 15 / exp * + (* π sqrt t 5 sqrt * 10 / erf * + (* 1 exp 20 nthroot t 2^^*

Numerical solutions

xoT xo transpose :

Method a: Trapezoidal rule

T Xa Xa' 13 mat D to te 12 mat xo 0.5 N DSolRes MCols :

Method b: Backward Euler

Tb Xb Xb' 13 mat D to te 12 mat xo 10 - N DSolRes MCols :

Method c: Forward Euler

Tc Xc Xc' 13 mat D to te 12 mat xo 0 N DSolRes MCols :

Method d: rkfixed from uni

Td Xd Xd' 13 mat xoT to te N D rkfixed MCols :

Method e: Rkadapt from uni

Te Xe Xe' 13 mat xoT to te N D Rkadapt MCols :

t x'.s maple x.s t diff feval : 1 exp 20 nthroot t 2^^t 5 sqrt * 25 sqrt 1 t 4 t + (* 20 / - exp - (* π sqrt 5 sqrt 2 t + (* 10 / erf - 15 sqrt / erf + (* 15 / exp * + (* π sqrt t 5 sqrt * 10 / erf * + (* 102 t + (* t 4 t + (* 20 / - exp * + 2015 / exp * 2 t + (2^20 / - exp * - 10 t 2^20 / - exp * + (* 10 /

(T's are the same)

T Td - normi T Te - normi 12 mat normi 0

Compare with the symbolic solution

k 1 N 1 + (range :

As exècted, Rkadapt it's better.

Xs k el T k el x.s :

Xs' k el T k el h + x.s Xs k el - h / :

Xs' k el T k el x'.s :

Xa Xs - normi 1.42 102 -^*

Xa' Xs' - normi 3.38 102 -^*

Xb Xs - normi 5.33 101 -^*

Xb' Xs' - normi 2.22 101 -^*

Xc Xs - normi 6.4 101 -^*

Xc' Xs' - normi 2.96 101 -^*

Xd Xs - normi 6.19 102 -^*

Xd' Xs' - normi 4.04 102 -^*

Xe Xs - normi 1.67 1010 -^*

Xe' Xs' - normi 2.66 102 -^*

T Xs augment T Xa o 3 black augment T Xb o 3 red augment T Xc o 3 green augment T Xd o 3 blue augment 5 1 sys T Xs' augment T Xa' o 3 red augment T Xb' o 3 black augment T Xc' o 3 green augment T Xd' o 3 blue augment 5 1 sys <content> DSolRes Example </content>

Example: How to build an example

Obviously, start from the solution:

Thus

n 14 range :

t S t t 2^sin * 2 + 1 t + ln 4 ln - e t -^t 2^cos 41 mat :

t x D dx n el t S n el t diff : (: 2 t 2^* t 2^cos * t 2^sin + 11 t + / 1 e t^/ - 2 t * t 2^sin * - 41 mat

You can use this D for, but just complicate the things

t x D 2 t 2^* x 4 el * x 1 el 2 - x 3 el ln / - 14 e x 2 el^* / x 3 el - 2 x 1 el 2 - (* - 41 mat :

Try to disable this for check that both D's are the same.

Given

to 0 :

te 2 :

N 25 :

Avoid t = 0

Now initials conditios are just

xo to S transpose : 2 1.3863 - 1114 mat

xoT xo transpose :

Numerical solutions

Method a: Trapezoidal rule

T Xa Ya Za Ua 15 mat D to te 12 mat xo 0.5 N DSolRes MCols :

Method b: Backward Euler

Tb Xb Yb Zb Ub 15 mat D to te 12 mat xo 10 - N DSolRes MCols :

Method c: Forward Euler

Tc Xc Yc Zc Uc 15 mat D to te 12 mat xo 0 N DSolRes MCols :

Method d: rkfixed from uni

Td Xd Yd Zd Ud 15 mat xoT to te N D rkfixed MCols :

Method e: Rkadapt from uni

Te Xa Ya Za Ua 15 mat xoT to te N D Rkadapt MCols :

Compare with the (given) symbolic solution S(t)

Hard to decide wich's better.

t S 1 el T Xa o 3 red augment T Xb o 3 black augment T Xc o 3 green augment T Xd o 3 blue augment 5 1 sys t S 2 el T Ya o 3 red augment T Yb o 3 black augment T Yc o 3 green augment T Yd o 3 blue augment 5 1 sys t S 3 el T Za o 3 red augment T Zb o 3 black augment T Zc o 3 green augment T Zd o 3 blue augment 5 1 sys t S 4 el T Ua o 3 red augment T Ub o 3 black augment T Uc o 3 green augment T Ud o 3 blue augment 5 1 sys

Alvaro