10.3.7 Approximating solutions of nonlinear second-order boundary value problems: bvpsolve
The bvpsolve command finds an approximate solution of a boundary value
problem
y″=f(x,y,y′), y(a)=α, y(b)=β
|
on an interval [a,b]. The procedure uses the method of nonlinear
shooting which is based on Newton and Runge-Kutta methods. Values of
y and its first derivative y′ are approximated at points
xk=a+k δ, where δ=b−a/N and k=0,1,…,N.
For the numeric tolerance (precision) threshold, the algorithm uses
epsilon specified in the session settings in Xcas (see
secrefssec:confcomp, item 9).
-
bvpsolve takes three mandatory arguments and four
optional arguments:
-
f(x,y,y′), an expression involving an independent variable
x, a dependent variable y and y′.
- [x=a..b,y], a list specifying the independent variable x,
its range [a,b] and the sought function y
- [α,β], a list of the boundary values of y.
- Optionally, A, an initial guess for y′(a) (by default,
(β−α)/(b−a)).
- Optionally, N, an integer greater than or equal to 2 (by
default 100).
- output=type or
Output=type, the type of the output, where
type can be one of:
-
list
- diff
- piecewise
- spline
(By default list).
- limit=M, a positive integer for a limit for the
number of iterations before the procedure is stopped
(by default there is no limit).
- bvpsolve(f(x,y,y′),[x=a..b,y],[α,β] ⟨,A,N,
output=type,limit=M⟩)
returns, for the different output types:
-
list, a list of pairs [xk,yk] where yk≈ y(xk),
- diff, a list of lists [xk,yk,y′k], where y′k≈ y′(xk),
- piecewise, a piecewise linear interpolation of the
points (xk,yk).
- spline, a piecewise spline interpolation of the
points (xk,yk), based on the values y′k computed in the
process.
Note that the shooting method is sensitive to roundoff errors and may
fail to converge in some cases, especially when y is a rapidly
increasing function. In the absence of convergence or if the maximum
number of iterations is exceeded, bvpsolve returns
undef. However, if the output type is list or
piecewise and if N>2, a slower but more stable
finite-difference method (which approximates only the function y) is
tried first.
Sometimes setting an initial guess A for y′(a) to a suitable value
may help the shooting algorithm to converge or to converge faster.
Examples.
-
Solve the problem
y″= | | (32+2 x3−y y′), 1≤ x≤ 3
|
with
boundary conditions y(1)=17 and y(3)=43/3.
Use N=20, which gives an x-step of 0.01.
Input:
bvpsolve((32+2x^3-y*y’)/8,[x=1..3,y],[17,43/3],20)
The output is shown in Table 10.1 (the middle two columns)
alongside with the values y(xk) of the exact solution y=x2+16/x
(the fourth column).
k | xk | yk | y(xk) |
0 | 1.0 | 17.0 | 17.0 |
1 | 1.1 | 15.7554961579 | 15.7554545455 |
2 | 1.2 | 14.7733911821 | 14.7733333333 |
3 | 1.3 | 13.9977543159 | 13.9976923077 |
4 | 1.4 | 13.388631813 | 13.3885714286 |
5 | 1.5 | 12.9167227424 | 12.9166666667 |
6 | 1.6 | 12.5600506483 | 12.56 |
7 | 1.7 | 12.3018096101 | 12.3017647059 |
8 | 1.8 | 12.1289281414 | 12.1288888889 |
9 | 1.9 | 12.0310865274 | 12.0310526316 |
10 | 2.0 | 12.0000289268 | 12.0 |
11 | 2.1 | 12.0290719981 | 12.029047619 |
12 | 2.2 | 12.1127475278 | 12.1127272727 |
13 | 2.3 | 12.2465382803 | 12.2465217391 |
14 | 2.4 | 12.4266798825 | 12.4266666667 |
15 | 2.5 | 12.650010254 | 12.65 |
16 | 2.6 | 12.9138537834 | 12.9138461538 |
17 | 2.7 | 13.2159312426 | 13.2159259259 |
18 | 2.8 | 13.5542890043 | 13.5542857143 |
19 | 2.9 | 13.9272429048 | 13.9272413793 |
20 | 3.0 | 14.3333333333 | 14.3333333333 |
Table 10.1: approximate and true values of the function y=x2+16/x on [1,3] |
- Solve
with the boundary conditions y(1)=0 and y(2)=ln256.
Obtain the solution as a piecewise spline interpolation for N=10 and
estimate the absolute error err of the approximation using
the exact solution y=x3 lnx and the romberg command for
numerical integration. You need to explicitly set an initial guess
A for the value y′(1) because the algorithm fails to converge with the
default guess A=ln256≈ 5.545. Therefore let A=1
instead.
Input:
f:=(x^2*diff(y(x),x)^2-9*y(x)^2+4*x^6)/x^5:;vars:=[x=1..2,y]:; yinit:=[0,ln(256),1]:; p:=bvpsolve(f,vars,yinit,10,output=spline):; err:=sqrt(romberg((p-x^3*ln(x))^2,x=1..2))
Output:
Note that, if the output type was set to list or
piecewise, the solution would have been found even without
specifying an initial guess for y′(1) because the algorithm would
automatically apply the alternative finite-difference method, which
converges.