Previous Up Next

21.4.3  Approximating solutions of y′′=f(t,y,y′) with boundary conditions

The bvpsolve command finds an approximate solution of a boundary value problem

  y′′=f(t,y,y′),   y(a)=α, y(b)=β

on the 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 tk=a+k δ, where δ=ba/N and k=0,1,…,N. For the numeric tolerance (precision) threshold, the algorithm uses epsilon specified in the session settings in Xcas (see Section 2.5.7, item 2.5.7).

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 y′′=1/8 (32+2 t3y y′) with 1≤ t≤ 3 and boundary conditions y(1)=17 and y(3)=43/3. Use N=20, which gives an t-step of 0.01.

bvpsolve((32+2t^3-y*y')/8,[t=1..3,y],[17,43/3],20)

The output is shown in Table 21.1 (the middle two columns) alongside with the values y(tk) of the exact solution y(t)=t2+16/t (the fourth column).


ktkyky(tk)
01.017.017.0
11.115.755496157915.7554545455
21.214.773391182114.7733333333
31.313.997754315913.9976923077
41.413.38863181313.3885714286
51.512.916722742412.9166666667
61.612.560050648312.56
71.712.301809610112.3017647059
81.812.128928141412.1288888889
91.912.031086527412.0310526316
102.012.000028926812.0
112.112.029071998112.029047619
122.212.112747527812.1127272727
132.312.246538280312.2465217391
142.412.426679882512.4266666667
152.512.65001025412.65
162.612.913853783412.9138461538
172.713.215931242613.2159259259
182.813.554289004313.5542857143
192.913.927242904813.9272413793
203.014.333333333314.3333333333
Table 21.1: Predicted and actual values of the function y(t)=t2+16/t on [1,3]

Solve y′′=t2 y2−9 y2+4 t6/t5 with 1≤ t≤ 2 and 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=t3 lnt 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.

f:=(t^2*diff(y(t),t)^2-9*y(t)^2+4*t^6)/t^5:; p:=bvpsolve(f,[t=1..2,y],[0,ln(256),1],10,output=spline):; err:=sqrt(romberg((p-t^3*ln(t))^2,t=1..2))
     
3.27751904973×10−6           

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.


Previous Up Next