<content> Utilities </content>

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

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

f.\003C\NFeval x.\003C\NFeval feval f.\003C\NFeval num2str (x.\003C\NFeval num2str) concat str2num 11 line :

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

S.NR Simple and fast Newton Raphson. fs a S.NR x a eval : ε ¦ε a UnitsOf * eval : 12 mat k 1 ¦MaxIters range y fs x feval : t x y ε * fs x ε + feval y - / - eval : 12 mat t x - abs ε < break x t : if 21 line for k ¦MaxIters ≤ t MaxIters rached error if 31 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 :

¦MaxIters 20 :

¦ε 107 -^:

X Y x CInterp UX X 1 el X 2 el - UnitsOf : UY Y 1 el Y 2 el - UnitsOf : 12 mat UY X UX / Y UY / x x UnitsOf / cinterp * x x UnitsOf / X UX / min X UX / max lele eval * 21 line :

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) = c2using 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) = c2using 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'NewtonMaxIters  rached. error if 81 line :

NewtonMaxIters 50 :

Boundary Values Problems for Beams

lbvp.2

x f x αβ N lbvp.2

solves the linear ODE

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

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

x a b 12 mat ≡

f is such that

x f x p x q x r 13 mat ≡

and

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

bvp.2

x y y' φ x Yo αβ 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.

From the free body diagram, we try to write

Short hands for plots X Y Y' xm Plot x m * y x m * 0 L lele (* x m * y' x m * 0 L lele (* X m / Y . 6 blue augment X m / Y' . 6 red augment xm xm y xm y' 21 mat o augment 51 sys 11 line :

y'' M E I * / ≡

Defaults

N ε 12 mat 50109 -^12 mat :

<content> Example beam </content>

Example

http://web.ncyu.edu.tw/~lanjc/lesson/C3/class/Chap06-A.pdf

Cantilever beam with uniformly distributed load

wo 400 lbf ft / * :

L 8 ft * :

E 29106^* psi * :

I 285 in 4^* :

(W12×35 shape)

FBD

x y'' wo 2 E * I * / x 2^* - :

Using RK Method

x y D y 2 el x y'' 21 mat :

IC:

y.0 0 :

y'.0 0 :

Sterilize units m kg s 13 mat 11113 mat :

X Y Y' 13 mat y.0 y'.0 stack L 0 N 1 - D rkfixed Cols :

Restore units m kg s Clear 1

X Y Y' 13 mat X m * Y m * Y' 13 mat :

Exact solution

x y wo 24 E * I * / x 4^- 4 L 3^* x * + 3 L 4^* - (* :

x y' x y x diff :

Maximun deflection

X Y Y' 0 Plot

Y N el in 0.0417 -

0 y in 0.0428 -

<content> Example beam </content>

Example

http://web.ncyu.edu.tw/~lanjc/lesson/C3/class/Chap06-A.pdf

Simple supported beam with distributed load distributed load of maximum intensity wo

wo 400 lbf ft / * :

L 8 ft * :

E 29106^* psi * :

I 285 in 4^* :

(W12×35 shape)

R.A wo L 4 / * :

R.C R.A :

FBD

x y'' wo 12 E * I * / 3 L 2^* x * 4 x 3^* - L / * :

Sterilize units m kg s 13 mat 11113 mat :

Using bvp.2

x y y' φ x y'' :

0 y 0 ≡ L y 0 ≡ 21 sys

αβ 10010023 mat :

X Y Y' 13 mat x y y' φ 0 L 12 mat 0 αβ N 1 - ε bvp.2 Cols :

or (better):

Using lbvp.2

x f 00 x y'' 13 mat :

X Y Y' 13 mat x f 0 L 12 mat αβ N 1 - lbvp.2 Cols :

Restore units m kg s Clear 1

X Y Y' 13 mat X m * Y m * Y' 13 mat :

Exact solution

x y wo x * 960 L * E * I * / 25 L 4^* 40 L 2^* x 2^* - 16 x 4^* + (* - :

x y' x y x diff :

Maximun deflection

X Y Y' 0.5 L * Plot

X Y 0.5 L * CInterp mm 0.068 -

0.5 L * y mm 0.0725 -

Only an approximate result.

<content> Example beam </content>

Example

http://web.ncyu.edu.tw/~lanjc/lesson/C3/class/Chap06-A.pdf

E 12 GPa * :

Q 300 N * :

x.B 2 m * :

L 3 m * :

I 40 mm * 80 mm * (3^* 12 / :

x y'' 13 / Q E I * / * x * x x.B ≤ 13 / Q E I * / * x * Q E I * / x x.B - (* - cases :

Sterilize units m kg s 13 mat 11113 mat :

Using bvp.2

x y y' φ x y'' :

0 y 0 ≡ L y 0 ≡ 21 sys

αβ 10010023 mat :

X Y Y' 13 mat x y y' φ 0 L 12 mat 0 αβ N 1 - ε bvp.2 Cols :

Using shooting method

x y D y 2 el x y'' 21 mat :

IC:

y.0 0 :

Guess y'.0 0.1 - :

BC:

y.L 0 :

y'.0 sol y.0 y'.0 stack 0 L N 1 - D rkfixed :

y'.0 Eq y'.0 sol N 2 el y.L - :

y'.0 Eq y'.0 S.NR :

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

Restore units m kg s Clear 1

Ξ Ψ Ψ' 13 mat Ξ m * Ψ m * Ψ' 13 mat :

x y Ξ Ψ x CInterp :

x y' Ξ Ψ' x CInterp :

Maximun deflection

x.max y' x.B S.NR : 1.621 m *

X Y Y' x.max Plot

x.max y mm 7.1271 -

<content> Example </content>

Example

https://edu.epito.bme.hu/local/coursepublicity/mod/resource/view.php?id=59032

L 4 m * :

EI 1.4 107^* N m 2^/ * :

q 10103^* N m / * :

x y'' 12 EI * / 1 y' 2^+ (3^sqrt * q * L x * x 2^- (* :

Sterilize units m kg s 13 mat 11113 mat :

Using bvp.2

x y y' φ x y'' :

0 y 0 ≡ L y' 0 ≡ 21 sys

M 10010023 mat :

X Y Y' 13 mat x y y' φ 0 L 12 mat 0 M N 1 - ε bvp.2 Cols :

Using shooting method

x y D y 2 el x y 1 el y 2 el φ 21 mat :

y.0 y.L 12 mat 0012 mat :

Guess y'.0 1 :

y'.0 sol y.0 y'.0 21 mat 0 L N 1 - D rkfixed :

y'.0 Eq y'.0 sol N 2 el y.L - :

y'.0 Eq y'.0 S.NR :

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

Restore units m kg s Clear 1

Ξ Ψ Ψ' 13 mat Ξ m * Ψ m * Ψ' 13 mat :

x y Ξ Ψ x CInterp :

x y' Ξ Ψ' x CInterp :

Maximun deflection

X Y Y' 0.5 L * Plot

x.max y' 0.4 L * S.NR : 1.9864 m *

x.max y mm 2.3795 -

Restore units m kg s Clear 1

Alvaro