Utilities

e5b2d458-1a86-4751-ab75-1a5708c59e5f 149 4 5 false false 0 decimal

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

Utilities

M Cols C 0 : k 1 M cols range : 12 mat C 1 k el M k col : 21 line :

A D C B TriDiagSys n B length : v n 1 matrix : y v : 13 mat w A 1 el 1 -^: y 1 el B 1 el w * : 12 mat k 2 n range v k 1 - el C k 1 - el w * eval : w A k el D k el v k 1 - el * - (1 -^eval : y k el B k el D k el y k 1 - el * - (w * eval : 31 line for k n 1 - (1 range y k el y k el v k el y k 1 + el * - eval : for y 51 line :

x y stem S 03 matrix : L 02 matrix : 12 mat k 1 x length range S S x k el y k el o 13 mat stack : L L x k el 0 x k el y k el x k el 032 mat stack : 21 line for L eval 31 line :

bvp.2

lbvp.2(f(x),[a,b],M,n) solves the linear BVP y'' + p(x)y' + q(x)y = φ(x) in a < x < b with
    f(x) = [p(x), q(x), φ(x)],  M = [α1 β1 c1; α2 β2 c2]  and 
    α1*y(a) + β1*y'(b) = c1   and  α2*y(a) + β2*y'(b) = c2
using central differences. 1 f x M n lbvp.2 h x 2 el x 1 el - n / : X x 1 el 0 n range h * + : ν n 1 + : P 0 : Q 0 : Φ 0 : α1 β1 c1 α2 β2 c2 23 mat M : 17 mat r 1 ν range P r el Q r el Φ r el 13 mat X r el f eval : for A D C B 14 mat 2 - h 2^Q * + (2 h P * - 2 / (2 h P * + 2 / (h 2^Φ * (14 mat eval : A 1 el D 1 el C 1 el B 1 el 14 mat β1 0 ≡ 100 c1 α1 / 14 mat A 1 el 2 h * α1 β1 / * D 1 el * + (02 B 1 el 2 h * c1 β1 / * D 1 el * + (14 mat if eval : A ν el D ν el C ν el B ν el 14 mat β2 0 ≡ 100 c2 α2 / 14 mat A ν el 2 h * α2 β2 / * C ν el * - (20 B ν el 2 h * c2 β2 / * C ν el * - (14 mat if eval : Y A D C B TriDiagSys : Y' Y : j 2 n range : 13 mat Y' j el Y j 1 + el Y j 1 - el - 2 h * / : Y' 1 el β1 0 ≡ Y 2 el Y 1 el - h / c1 β1 / α1 β1 / Y 1 el * - if : Y' ν el β2 0 ≡ Y ν el Y n el - h / c2 β2 / α2 β2 / Y ν el * - if : 13 mat X Y Y' augment eval 81 line :

bvp.2(f(x),[a,b],M,n) solves the linear BVP y'' = f(x,y,y') in a < x < b with
    f(x) = [p(x), q(x), φ(x)] , M = [α1 β1 c1; α2 β2 c2]  and 
    α1*y(a) + β1*y'(b) = c1   and  α2*y(a) + β2*y'(b) = c2
using central differences. 3 f x Y M n ε bvp.2 h x 2 el x 1 el - n / : X x 1 el 0 n range h * + : ν n 1 + : P 0 : Q 0 : Φ 0 : α1 β1 c1 α2 β2 c2 23 mat M : 17 mat Y 0 ≡ ψ x 1 el α1 * β1 + α1 x 2 el α2 * β2 + α2 22 mat (1 -^c1 c1 21 mat * : Y ψ 2 el ψ 1 el X * + : 12 mat 0 if Δy Y max Y min - ν / (2^h 2^+ sqrt eval : Δy' Δy h 1 -^* : 12 mat x y y' fy x y Δy + y' f x y Δy - y' f - 2 Δy * / : x y y' fy' x y y' Δy' + f x y y' Δy' - f - 2 Δy' * / : A ν 1 matrix : D A : C A : B A : Y' A : j 2 n range : 16 mat Iter 1 NewtonMaxIters range Y' j el Y j 1 + el Y j 1 - el - 2 h * / : Y' 1 el β1 0 ≡ Y 2 el Y 1 el - h / c1 β1 / α1 β1 / Y 1 el * - if : Y' ν el β2 0 ≡ Y ν el Y n el - h / c2 β2 / α2 β2 / Y ν el * - if : 13 mat r 1 ν range P r el Q r el Φ r el 13 mat X r el Y r el Y' r el fy X r el Y r el Y' r el fy' X r el Y r el Y' r el f 13 mat eval : for A D C 13 mat 2 - h 2^P * - (2 h Q * - 2 / (2 h Q * + 2 / (13 mat eval : B j el 2 Y j el * - Y j 1 + el + Y j 1 - el + h 2^Φ j el * - eval : A 1 el D 1 el C 1 el B 1 el 14 mat β1 0 ≡ 100 Y 1 el c1 α1 / - 14 mat A 1 el 2 h * α1 β1 / * B 1 el * + (022 Y 2 el * 21 - h α1 β1 / * + (* Y 1 el * + 2 h * c1 β1 / * - h 2^Φ 1 el * - (14 mat if eval : A ν el D ν el C ν el B ν el 14 mat β2 0 ≡ 100 Y ν el c2 α2 / - 14 mat A ν el 2 h * α2 β2 / * C ν el * - (202 Y n el * 21 h α2 β2 / * + (* Y ν el * - 2 h * c2 β2 / * + h 2^Φ ν el * - (14 mat if eval : Yo A D C B - TriDiagSys : Y Y Yo + eval : 12 mat Yo norme ε Y norme * < break continue if 81 line for Iter NewtonMaxIters < X Y Y' augment eval \0027\NewtonMaxIters  rached. error if 81 line :

NewtonMaxIters 50 :

Boundary Values Problem

lbvp.2

x f x M N lbvp.2

solves the linear ODE

y'' p y' * + q y * + r ≡

x a b 12 mat ≡

subject to the boundary conditions

α 1 el a y * β 1 el a y' * + c 1 el ≡ α 2 el b y * β 2 el b y' * + c 2 el ≡ 21 sys

f is such that

x f x p x q x r 13 mat ≡

and

M α1 β1 c1 α2 β2 c2 23 mat ≡

bvp.2

x y y' φ x 0 M N ε bvp.2

solves the non linear ODE

y'' x y y' φ ≡

x a b 12 mat ≡

subject to the the same boundary conditions, with Yo as guess for the solution with dimension N+1.. If Yo≡0, bvp.2 try with a line between f(a) and f(b) as guess.. ε is used as the tolerance for a Newton solver.

Short hands for plots

Short hands for plots

x xb X Y Y' x y x y' G Ys Y's Clear 1

Err N X length : n 1 N range : Ys n el X n el y : 13 mat Y Ys - norme 21 line :

PlotXY N X length : x y X Y . 8 augment a b 12 mat Y 1 el Y N el 12 mat stem a b 12 mat transpose Y 1 el Y N el 12 mat transpose o augment 41 sys 21 line :

PlotYY' N X length : n 1 N range : 12 mat Ys n el X n el y : Ys' n el X n el y' : augment Y Y' . 8 augment Y 1 el Y N el 12 mat Y' 1 el Y' N el 12 mat stem Y 1 el Y N el 12 mat transpose Y' 1 el Y' N el 12 mat transpose o augment 41 sys 21 line :

lbvp example

Example

Solve

y'' 2 h * y' * + h 2^h' + (y * + 0 ≡

a b 12 mat 0512 mat :

http://eqworld.ipmnet.ru/en/solutions/ode/ode0244.pdf

subject to the boundary conditions

bc a y 3 a y' * + 0.1 ≡ b y' 6 - ≡ 21 sys :

Numeric solution

h x x cos * :

h' h x diff :

x f 2 h * h 2^h' + 013 mat :

M 13 0.1 016 - 23 mat :

Analytic solution

x y A B x * + (e x x sin * x cos + (-^* :

x y' x y x diff :

N ε 12 mat 50103 -^12 mat :

A B 21 mat bc sys2mat.1 A B 21 mat roots :

Results

X Y Y' 13 mat x f a b 12 mat M 100 lbvp.2 Cols :

Err 0.16

PlotXY PlotYY'

Note

Another way to solve for the integration constants is

A B Clear 1

eqM M 11 el a y * M 12 el a y' * + M 13 el - M 21 el b y * M 22 el b y' * + M 23 el - 21 mat :

eqM A B 21 mat roots 0.7939 0.174 - 21 mat

x y x y' A B Clear 1

lbvp example

Example

Solve

y'' 3 y' * + 4 y * - 0 ≡

a b 12 mat 0112 mat :

subject to the boundary conditions

0 y 1 ≡ 1 y 1 ≡ 21 sys

Numeric solution

x f 34 - 013 mat :

M 10110123 mat :

Analytic solution

x y e 44 x * - (^e x^+ e x 1 + (^+ e x 2 + (^+ e x 3 + (^+ (1 e + e 2^+ e 3^+ e 4^+ (/ :

x y' x y x diff :

Results

X Y Y' 13 mat x f a b 12 mat M 50 lbvp.2 Cols :

Err 0

PlotXY PlotYY'

x x y x y' Clear 1

lbvp example

Example

Solve

y'' 3 y' * + 6 + 5 ≡

a b 12 mat 1312 mat :

subject to the boundary conditions

1 y 3 ≡ 3 y 23 y' * + 5 ≡ 21 sys

Numeric solution

x f 36513 mat :

M 10312523 mat :

Analytic solution

x y 0.0696737 e 1.5 x * -^* 11.9605 e 1.5 x *^* 1243.5 1.93649 x * sin * + 2857.67 1.93649 x * cos * + (* :

x y' x y x diff :

Results

X Y Y' 13 mat x f a b 12 mat M 100 lbvp.2 Cols :

Err 0.21

PlotXY PlotYY'

x y x y' Clear 1

lbvp example

Example

Solve

y'' 3 p * p x 2^+ (2^/ y * + 0 ≡

a b 12 mat 0.1 - 0.1 12 mat :

subject to the boundary conditions

a y y1 ≡ b y y2 ≡ 21 sys

https://viewer.mathworks.com/?viewer=plain_code&url=https%3A%2F%2Fwww.mathworks.com%2Fmatlabcentral%2Fmlc-downloads%2Fdownloads%2Fsubmissions%2F3819%2Fversions%2F6%2Fcontents%2FBVP_tutorial%2FBVP_examples%2Fex2bvp.m&embed=web

Analytic solution

p 105 -^:

x y x p x 2^+ sqrt / :

x y' x y x diff :

Numeric solution

x f 03 p * p x 2^+ (2^/ 013 mat :

M 10 a y 10 b y 23 mat :

Results

X Y Y' 13 mat x f a b 12 mat M 500 lbvp.2 Cols :

< Big N

Err 0.29

PlotXY PlotYY'

x y x y' p Clear 1

lbvp example

Example

Solve

y'' 2 x / y' * + y x 4^/ + 0 ≡

a b 12 mat 13 π * / 3 π / 12 mat :

subject to the boundary conditions

13 π * / y 0 ≡ 3 π / y 3 sqrt 2 / ≡ 21 sys

https://www.mathworks.com/help/matlab/ref/bvp4c.html

Numeric solution

x f 2 x / 1 x 4^/ 013 mat :

M 100103 sqrt 2 / 23 mat :

Analytic solution

x y 1 x / sin :

x y' x y x diff :

Results

X Y Y' 13 mat x f a b 12 mat M 200 lbvp.2 Cols :

Err 0.11

PlotXY PlotYY'

x y x y' Clear 1

bvp with two solutions

Example

Solve

y'' y abs + 0 ≡

a b 12 mat 0412 mat :

subject to the boundary conditions

a y 0 ≡ b y 2 - ≡ 21 sys

http://www.ece.northwestern.edu/local-apps/matlabhelp/techdoc/ref/bvp4c.html

This problem have two solutions. I compare the bvp solution with the shooting method, because can't found a symbolic expression.

Numeric solution

x y y' φ y abs - :

M 100102 - 23 mat :

X Y Y' 13 mat x y y' φ a b 12 mat 0 M N ε bvp.2 Cols :

Shoot solution

y.a y.b 12 mat Y 1 el Y N 1 + el 12 mat :

Guess:

y'.a Y' 1 el :

x y D y 2 el y 1 el abs - 21 mat :

y'.a sol y.a y'.a 21 mat a b N D al_rkckadapt :

y'.a Eq y'.a sol N 1 + 2 el y.b - :

y'.a y'.a Eq al_nleqsolve 1 el : 0.0733 -

Ξ Ψ Ψ' 13 mat y'.a sol Cols :

x y Ξ Ψ x cinterp :

x y' Ξ Ψ' x cinterp :

Results

Err 0

PlotXY PlotYY'

For the other solution: call bvp with a new guess

X Y Y' 13 mat x y y' φ a b 12 mat Yo 1 N 1 + (range el 1 : M N ε bvp.2 Cols :

Solve the shoot with a new guess from bvp

y'.a Y' 1 el :

y'.a y'.a Eq al_nleqsolve 1 el : 2.0666

Ξ Ψ Ψ' 13 mat y'.a sol Cols :

x y Ξ Ψ x cinterp :

x y' Ξ Ψ' x cinterp :

Results

Err 0.01

PlotXY PlotYY'

x y x y' Clear 1

Shooting method

Shooting Method

Solve this BVP using the Shooting method

ode x y x 2 diff 12 / x y x diff * + x y + 5 - :

a b 12 mat 1 - 612 mat :

y.a y.b 12 mat 7312 mat :

Convert to a system and solve with an ode solver

x y D y 2 el y 2 el 2 / - y 1 el - 5 + 21 mat :

N 100 :

y'.a sol y.a y'.a 21 mat a b N D al_rkckadapt :

Equation to solve

y'.a Eq y'.a sol N 1 + 2 el y.b - :

Guess for a solution

y'.a 10 - :

Solve

y'.a y'.a Eq al_nleqsolve 1 el : 27.5704 -

X Y Y' 13 mat y'.a sol Cols :

lbvp solution

x f 12 / 1513 mat :

M 10710323 mat :

Ξ Ψ Ψ' 13 mat x f a b 12 mat M N lbvp.2 Cols :

Results

Y Ψ - norme 0.34

Ξ Ψ augment X Y . augment y'.a x a - ( * y.a + 3 1 sys Ψ Ψ' augment Y Y' . augment y.a y'.a o augment 3 1 sys

Alvaro