bvp .

“msharpmath”, The Simple is the Best

[100] 007 Tutorial bvp            boundary value problems

Click the above button for full-text tutorial of  ‘bvp’.

//==================================================================================
//  Umbrella  ‘bvp’  for Boundary Value Problems
//==================================================================================

//————————————————————————————————————
//  Standard Differential Equation
//————————————————————————————————————
cz*{f””}
+cw*{f”’} + cfw*{ff”’} + cuw*{f’f”’} + cvw*{f”f”’} + cww*{f”’f”’}
+cv*{f”}  + cfv*{ff”}  + cuv*{f’f”}  + cvv*{f”f”}
+cu*{f’}   + cfu*{ff’}   + cuu*{f’f'}
+cf*{f}    + cff*{ff} = cs

//————————————————————————————————————
//  Standard Boundary Conditions
//————————————————————————————————————

bw*{f”’} + bv*{f”} + bu*{f’} + bf*{f} = bs

In an interval (a,b),

bvp .x(a,b) …

can have BCs the numbers of which are

one at x=a, one at x=b for 2nd-order ODE
two at x=a, one at x=b for 3rd-order ODE
two at x=a, two at x=b for 4th-order ODE

Other types of BCs cannot be handled yet. Especially for the 3rd-order BVP,
the starting point can be switched by the ending point to meet the above condition such that

bvp .x(b,a) …

so that an additional type (one at x=a, two at x=b) can be handled.

//————————————————————————————————————
//  Spokes
//————————————————————————————————————
.peep                                       // peep the iteration procedure
.relax(relax=0.5)                  // under-relaxation factor
.tol/abstol(abstol=1.e-5)    // absolute tolerance
.reltol(reltol=1.e-5)              // relative tolerance
.maxiter(maxiter=500)      // maximum iteration count
.shift(s)                                  // eigenvalue shift
.plot                                        // plot default data
.plot(f,f’,f”,…)                        // plot user-defined data
.return(f,f’,f”,…)                   // return user-defined data
.togo(Y=y,Y1=y’,…)             // takeout each field

//————————————————————————————————————
//  Examples
//————————————————————————————————————

%> simple BVP =======================================


#> bvp .x(0,pi) ( {y”}+{y} = 0, [ {y}=0, {y'}=-1 ] ).plot;

%> shock near at x = 0 =======================================


#> eps = 1.e-4;;
#{
bvp .x(-1,1) (       // bvp .x[101,1](-1,1) equivalently
eps*{y”}+x*{y’} = eps*pi^2*cos(pi*x)-pi*x*sin(pi*x), [ {y} = -2, {y} = 0 ]
)
.plot .return(x,y,y’);
#}

ans =
[            -1            -2    0.00395044 ]
[         -0.98      -1.99799      0.197356 ]
[         -0.96      -1.99203      0.398039 ]
…..
[          0.96    0.00796702     -0.398039 ]
[          0.98    0.00201307     -0.197356 ]
[             1  8.36201e-009   -0.00395044 ]


%> droplet shape ===========================================


#{
bvp .x[21](-1,1) ( {h”} + (1+h’^2)^1.5  – (1+h’^2)^1.5 *{h} = 0, [ {h} = 0, {h} = 0 ] )
.plot
.return(x,h,h’);
#}
ans =
[            -1  4.18101e-009       1.05936 ]
[          -0.9     0.0936755       0.81415 ]

[          -0.1      0.387533     0.0611795 ]
[             0      0.390592  1.72213e-013 ]
[           0.1      0.387533    -0.0611795 ]

[           0.9     0.0936755      -0.81415 ]
[             1  4.18101e-009      -1.05936 ]

%> Falkner-Skan equations ====================================


#{
.hold;
for.b(0,1, 0.1) {
bvp .x[101](0,8) ( {f”’}+{ff”}+b-b*{f’f'} = 0, [ {f} = 0, {f'} = 0, {f'} = 1 ] )
.plot(f’,f”)
.return(f”)(1);
}
plot;
#}

ans =      0.46983977
ans =      0.58729653
ans =      0.68700019
ans =      0.77508197
ans =      0.85478704
ans =      0.92808654
ans =      0.99628555
ans =       1.0603011
ans =       1.1208068
ans =       1.1783141
ans =       1.2332224

%> Coupled BVP ====================================
#{
bvp.x[21](0,8) { f,g } (
{f”’}+3*{ff”}-2*{f’f'} + g = 0,  [ {f} = 0, {f'} = 0, {f} = 0 ],
{g”}+3*f*{g’} = 0,                [ {g} = 1, {g} = 0 ]
) .plot ;
#}

//————————————————————————————————————
//  Eigenvalue Problems
//————————————————————————————————————
Standard forms are

{u”} + {%u} = 0     // lambda is denoted by %

Note that non-solitary {%u} such as 2*{%u} is not allowed

Eigenvalue shifting by Spoke .shift(s)

%> basic example
#> x = (0,pi).span(101);;

#> bvp[x]( {u”}+{%u} = 0, [ {u}=0, {u}=0 ] );
ans =      0.99991775

 

%> second-order example


#> lam =
bvp .x[101](0,pi) (
{y”} + {%y} = 0,
[ {y}-{y'} = 0, {y} = 0 ]
).plot ; // dependent variable y is assumed for print

#> lam;
lam =      0.62142265

//————————————————
//  shifting
//————————————————
Spoke .shift(s) is applied to find eigenvalue nearest to s

#> x = (0,pi).span(101);;

#> lam = bvp[x]( {y”}+{%y} = 0, [ {y}=0, {y}=0 ] ) .shift(5);
lam =       3.9986842

#> lam = bvp[x]( {y”}+{%y} = 0, [ {y}=0, {y}=0 ] ) .shift(10);
lam =         8.99334

#> lam = bvp.x[500](0,pi)( {y”}+{%y} = 0, [ {y}=0, {y}=0 ] ) .shift(10);
lam =       8.9997325

#> lam = bvp.x[1500](0,pi)( {y”}+{%y} = 0, [ {y}=0, {y}=0 ] ) .shift(10);
lam =       8.9999704

Exact values are lam = 1^2, 2^2, 3^2, 4^2, … , n^2. A dense grid enhances the accuracy of the eigenvalue found.

//————————————————
//  bending of a beam (4th-order ODE)
//————————————————


#> lam = bvp.x(0,8) (
-1000*{u””} + {%u} = 0,
[ {u}=0, {u'}=0, {u''}=0, {u'''}=0 ]
)
.plot(u,u’,u”,u”’);

#> sqrt(lam);
ans =       1.7374164

Comments are closed.