1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
|
bvode Scilab Group Scilab Function bvode
NAME
bvode - boundary value problems for ODE
CALLING SEQUENCE
[z]=bvode(points,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...
fsub1,dfsub1,gsub1,dgsub1,guess1)
PARAMETERS
z The solution of the ode evaluated on the mesh given by points
points an array which gives the points for which we want the solution
ncomp number of differential equations (ncomp <= 20)
m a vector of size ncomp. m(j) gives the order of the j-th
differential equation
( mstar = m(1) + ... + m(ncomp) <= 40 )
aleft left end of interval
aright right end of interval
zeta zeta(j) gives j-th side condition point (boundary point). must
have
zeta(j) <= zeta(j+1)
all side condition points must be mesh points in all meshes
used, see description of ipar(11) and fixpnt below.
ipar an integer array dimensioned at least 11. a list of the
parameters in ipar and their meaning follows some parameters
are renamed in bvode; these new names are given in parentheses.
ipar(1) ( = nonlin )
= 0 if the problem is linear
ipar(2) = number of collocation points per subinterval (= k
) where
max m(i) <= k <= 7 .
if ipar(2)=0 then bvode sets
k = max ( max m(i)+1, 5-max m(i) )
ipar(3) = number of subintervals in the initial mesh ( = n
). if ipar(3) = 0 then bvode arbitrarily sets n = 5.
ipar(4) = number of solution and derivative tolerances. ( =
ntol ) we require
0 < ntol <= mstar.
ipar(5) = dimension of fspace ( = ndimf ) a real work array.
its size provides a constraint on nmax. choose
ipar(5) according to the formula
ipar(5) >= nmax*nsizef
where
nsizef = 4 + 3 * mstar + (5+kd) * kdm +
(2*mstar-nrec) * 2*mstar.
ipar(6) = dimension of ispace ( = ndimi )an integer work
array. its size provides a constraint on nmax, the
maximum number of subintervals. choose ipar(6)
according to the formula
ipar(6) >= nmax*nsizei
where
nsizei = 3 + kdm
with
kdm = kd + mstar ; kd = k * ncomp ;
nrec = number of right end boundary conditions.
ipar(7) output control ( = iprint )
= -1 for full diagnostic printout
= 0 for selected printout
= 1 for no printout
ipar(8) ( = iread )
= 0 causes bvode to generate a uniform initial mesh.
= xx Other values are not implemented yet in Scilab
= 1 if the initial mesh is provided by the user.
it is defined in fspace as follows: the
mesh
aleft=x(1)<x(2)< ... <x(n)<x(n+1)=aright
will occupy fspace(1), ..., fspace(n+1). the
user needs to supply only the interior mesh
points fspace(j) = x(j), j = 2, ..., n.
= 2 if the initial mesh is supplied by the user
as with ipar(8)=1, and in addition no adaptive
mesh selection is to be done.
ipar(9) ( = iguess )
= 0 if no initial guess for the solution is
provided.
= 1 if an initial guess is provided by the user in
subroutine guess.
= 2 if an initial mesh and approximate solution
coefficients are provided by the user in fspace.
(the former and new mesh are the same).
= 3 if a former mesh and approximate solution
coefficients are provided by the user in fspace,
and the new mesh is to be taken twice as coarse;
i.e.,every second point from the former mesh.
= 4 if in addition to a former initial mesh and
approximate solution coefficients, a new mesh is
provided in fspace as well. (see description of
output for further details on iguess = 2, 3, and
4.)
= 0 if the problem is regular
= 1 if the first relax factor is =rstart, and the
nonlinear iteration does not rely on past
covergence (use for an extra sensitive nonlinear
problem only).
= 2 if we are to return immediately upon (a) two
successive nonconvergences, or (b) after
obtaining error estimate for the first time.
ipar(11) = number of fixed points in the mesh other than
aleft and aright. ( = nfxpnt , the dimension of
fixpnt) the code requires that all side condition
points other than aleft and aright (see description
of zeta ) be included as fixed points in fixpnt.
ltol an array of dimension ipar(4). ltol(j) = l specifies that the
j-th tolerance in tol controls the error in the l-th
component of z(u). also require that
1<=ltol(1)<ltol(2)< ... <ltol(ntol)<=mstar
tol an array of dimension ipar(4). tol(j) is the error tolerance
on the ltol(j) -th component of z(u). thus, the code attempts
to satisfy for j=1,...,ntol on each subinterval
abs(z(v)-z(u)) \le tol(j)*abs(z(u)) +tol(j)
ltol(j) ltol(j)
if v(x) is the approximate solution vector.
fixpnt an array of dimension ipar(11). it contains the points, other
than aleft and aright, which are to be included in every mesh.
externals The function fsub,dfsub,gsub,dgsub,guess are Scilab externals
i.e. functions (see syntax below) or the name of a Fortran
subroutine (character string) with specified calling sequence
or a list. An external as a character string refers to the
name of a Fortran subroutine. The Fortran coded function
interface to bvode are specified in the file fcol.f.
fsub name of subroutine for evaluating
t
f(x,z(u(x))) = (f ,...,f )
1 ncomp
at a point x in (aleft,aright). it should have the heading
[f]=fsub(x,z) where f is the vector containing the value
of fi(x,z(u)) in the i-th component and
t
z(u(x))=(z(1),...,z(mstar))
is defined as above under purpose .
dfsub
name of subroutine for evaluating the Jacobian of f(x,z(u)) at
a point x. it should have the heading [df]=dfsub (x , z )
where z(u(x)) is defined as for fsub and the (ncomp) by
(mstar) array df should be filled by the partial
derivatives of f, viz, for a particular call one
calculates
df(i,j) = dfi / dzj, i=1,...,ncomp
j=1,...,mstar.
gsub name of subroutine for evaluating the i-th component of
g(x,z(u(x))) = g (zeta(i),z(u(zeta(i))))
i
at a point x = zeta(i) where
1<=i<=mstar.
it should have the heading[g]=gsub (i , z) where z(u) is as
for fsub, and i and g=gi are as above. note that in
contrast to f in fsub , here only one value per call is
returned in g.
dgsub
name of subroutine for evaluating the i-th row of the Jacobian
of g(x,u(x)). it should have the heading [dg]=dgsub (i , z
) where z(u) is as for fsub, i as for gsub and the
mstar-vector dg should be filled with the partial
derivatives of g, viz, for a particular call one calculates
dg(i,j) = dgi / dzj, j=1,...,mstar.
guess
name of subroutine to evaluate the initial approximation for
z(u(x)) and for dmval(u(x))= vector of the mj-th
derivatives of u(x). it should have the heading [z,dmval]=
guess (x ) note that this subroutine is used only if
ipar(9) = 1, and then all mstar components of z and
ncomp components of dmval should be specified for any x,
aleft <= x <= aright .
DESCRIPTION
this package solves a multi-point boundary value problem for a mixed
order system of ode-s given by
(m(i))
u = f ( x; z(u(x)) ) i = 1, ... ,ncomp
i i
aleft < x < aright,
g ( zeta(j); z(u(zeta(j))) ) = 0 j = 1, ... ,mstar
j
mstar = m(1)+m(2)+...+m(ncomp),
where
t
u = (u , u , ... ,u ) is the exact solution vector
1 2 ncomp
(mi)
u is the mi=m(i) th derivative of u
i i
(1) (m1-1) (mncomp-1)
z(u(x)) = ( u (x),u (x),...,u (x),...,u (x) )
1 1 1 ncomp
f (x,z(u)) is a (generally) nonlinear function of
i
z(u)=z(u(x)).
g (zeta(j);z(u)) is a (generally) nonlinear function
j
used to represent a boundary condition.
the boundary points satisfy
aleft <= zeta(1) <= .. <= zeta(mstar) <= aright
the orders mi of the differential equations satisfy
1<=m(i)<=4.
EXAMPLE
deff('df=dfsub(x,z)','df=[0,0,-6/x**2,-6/x]')
deff('f=fsub(x,z)','f=(1 -6*x**2*z(4)-6*x*z(3))/x**3')
deff('g=gsub(i,z)','g=[z(1),z(3),z(1),z(3)];g=g(i)')
deff('dg=dgsub(i,z)',['dg=[1,0,0,0;0,0,1,0;1,0,0,0;0,0,1,0]';
'dg=dg(i,:)'])
deff('[z,mpar]=guess(x)','z=0;mpar=0')// unused here
deff('u=trusol(x)',[ //for testing purposes
'u=0*ones(4,1)';
'u(1) = 0.25*(10*log(2)-3)*(1-x) + 0.5 *( 1/x + (3+x)*log(x) - x)'
'u(2) = -0.25*(10*log(2)-3) + 0.5 *(-1/x^2 + (3+x)/x + log(x) - 1)'
'u(3) = 0.5*( 2/x^3 + 1/x - 3/x^2)'
'u(4) = 0.5*(-6/x^4 - 1/x/x + 6/x^3)'])
fixpnt=0;m=4;
ncomp=1;aleft=1;aright=2;
zeta=[1,1,2,2];
ipar=zeros(1,11);
ipar(3)=1;ipar(4)=2;ipar(5)=2000;ipar(6)=200;ipar(7)=1;
ltol=[1,3];tol=[1.e-11,1.e-11];
res=aleft:0.1:aright;
z=bvode(res,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...
fsub,dfsub,gsub,dgsub,guess)
z1=[];for x=res,z1=[z1,trusol(x)]; end;
z-z1
AUTHOR
u. ascher, department of computer science, university of british
columbia, vancouver, b. c., canada v6t 1w5 g. bader, institut f.
angewandte mathematik university of heidelberg im neuenheimer feld
294d-6900 heidelberg 1 Fotran subroutine colnew.f
SEE ALSO
fort, link, external, ode, dassl
|
