Scilab Function
Last update : 8/2/2006

bvodeS - simplified call of bvode

### Calling Sequence

  `z=bvodeS(x,m,n,a,b,fsub,gsub,zeta, [ystart,dfsub,dgsub,fixpnt,ndimf,ndimi,ltol,tol,ntol,nonlin, collpnt,subint,iprint,ireg,ifail])`

### Parameters

- `x`: array of points at which approximations of the solution will be computed. The points x(i) must be given in increasing order.
- `m`: array of orders for the differential equations given in `fsub`.
- `n`: number of differential equations, length of `m`.
- `a`: left end of solution interval, where a<=x(1).
- `b`: right end of solution interval, where x($)<=b.
- `fsub`: name of a Scilab function for evaluating the rhs of the differential equation system to be solved; fsub also may be a list for parameter transfer.
- `gsub`: name of a Scilab function for evaluating the boundary conditions for the given differential equation system; gsub may also be a list as for fsub.
- `zeta`: array of points at which the boundary or side conditions are given in increasing order. Each point must occur as often as boundary or side conditions are given there. Each side condition point other than a or b must be included in the array fixpnt. The length of zeta should be sum(m).
- `z`: array of all derivatives of the solution functions up to the order m(i)-1 for the i-th function. (See the examples below.) The length of z should be sum(m).
- `ystart,dfsub,dgsub,fixpnt,ndimf,ndimi,ltol,tol,ntol,nonlin, collpnt,subint,iprint,ireg,ifail`: These optional arguments may be called by name in any order in the form argument=name. The meaning of ystart to ireg is given in the `bvode` help page. The Scilab functions `dfsub` and `dgsub` for evaluating the Jacobians may also be lists for parameter transfer. The function `ystart` is called `guess` in bvode.
- `ifail`: if ifail=1, all parameters needed for the call of bvode are displayed.

### Description

This interface program simplifies the call of `bvode`, a program for the numerical solution of multi-point boundary value problems for mixed order systems of ordinary differential equations. The Scilab program bvode is adapted from the fortran program colnew. See the paper by U. Asher, J. Christiansen, and R.D. Russel: Collocation software for boundary-value ODE's, ACM Trans. Math. Soft. 7:209-222, 1981. The following examples should demonstrate not only how such systems can be solved with the help of bvodeS, but also should emphazise some important problems which occur with boundary value problems for ordinary ode's.

### Examples

    // 1. Modified example from help bvode.

    // DE:  y1''''(x)=(1-6*x*x*y1'''(x)-6*x*y1''(x))/(y2(x)^3)
    //         y2'(x)=1
    // z=[y1 ; y1' ; y1'' ; y1''' ; y2]
    // BV: y1(1)=0 y1''(1)=0
    // BV: y1(2)=1 y1''(2)=0 y2(2)=2

    function RhS=fsub(x,z)
      RhS=[(1-6*x*x*z(4)-6*x*z(3))/(z(5)^3);1]
    endfunction

    function g=gsub(i,z)
      g=[z(1) z(3) z(1)-1 z(3) z(5)-2]
      g=g(i)
    endfunction

    function [z,lhS]=ystart(x)
      z=zeros(5,1);z(5)=1;
      lhS=[0;1];
    endfunction

    n=2;
    m=[4 1];
    N=100;
    a=1; b=2;
    zeta=[a a b b b];
    x=linspace(a,b,N);

    ltol=4; // We want to change the default error for y1'''.
    tol=1e-12;
    tic()
    z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ltol=ltol,tol=tol,ystart=ystart);
    // Try tol=1e-14 etc.
    toc()

    function z=yex(x) // True solution
      z=zeros(5,1);
      z(1)=0.25*(10*log(2)-3)*(1-x)+0.5*(1/x+(3+x)*log(x)-x)+(x-1)
      z(2)=-0.25*(10*log(2)-3)+0.5*(-1/x^2+(3+x)/x+log(x)-1)+1
      z(3)=0.5*(2/x^3+1/x-3/x^2)
      z(4)=0.5*(-6/x^4-1/x/x+6/x^3)
      z(5)=x
    endfunction

    zex=[];for xx=x, zex=[zex yex(xx)]; end
    scf(0); clf();
    plot2d(x,abs(z-zex)',style=[1 2 3 5 6])
    xtitle('Absolute error','x',' ')
    legend(['z1(x)';'z2(x)';'z3(x)';'z4(x)';'z5(x)'])

    // example #2. An eigenvalue problem

    // y''(x)=-la*y(x)
    // BV: y(0)=y'(0); y(1)=0
    // Eigenfunctions and eigenvalues are y(x,n)=sin(s(n)*(1-x)), la(n)=s(n)^2,
    // where s(n) are the zeros of f(s,n)=s+atan(s)-(n+1)*pi, n=0,1,2,...
    // To get a third boundary condition, we choose y(0)=1
    // (With y(x) also c*y(x) is a solution for each constant c.)
    // We solve the following ode system:
    // y''=-la*y
    // la'=0
    // BV: y(0)=y'(0), y(0)=1; y(1)=0
    // z=[y(x) ; y'(x) ; la]

    function rhs=fsub(x,z)
      rhs=[-z(3)*z(1);0]
    endfunction

    function g=gsub(i,z)
      g=[z(1)-z(2) z(1)-1 z(1)]
      g=g(i)
    endfunction

    // The following start function is good for the first 8 eigenfunctions.
    function [z,lhs]=ystart(x,z,la0)
      z=[1;0;la0]
      lhs=[0;0]
    endfunction

    a=0;b=1;
    m=[2;1];
    n=2;
    zeta=[a a b];
    N=101;
    x=linspace(a,b,N)';

    // We have s(n)-(n+1/2)*pi -> 0 for n to infinity.
    la0=input('n-th eigenvalue: n= ?');la0=(%pi/2+la0*%pi)^2;

    z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,la0));

    clf()
    plot2d(x,[z(1,:)' z(2,:)'],style=[5 1],axesflag=5)
    xtitle(['Startvalue =  '+string(la0);'Eigenvalue = '+string(z(3,1))],'x',' ')
    legend(['y(x)';'y''(x)'])

    // example #3. A boundary value problem with more than one solution.

    // DE: y''(x)=-exp(y(x))
    // BV: y(0)=0; y(1)=0
    // This boundary value problem has more than one solution.
    // It is demonstrated how to find two of them with the help of
    // some preinformation of the solutions y(x) to build the function ystart.
    // z=[y(x);y'(x)]

    a=0;b=1;m=2;n=1;
    zeta=[a b];
    N=101;
    tol=1e-8*[1 1];
    x=linspace(a,b,N);

    function rhs=fsub(x,z),rhs=-exp(z(1));endfunction

    function g=gsub(i,z)
      g=[z(1) z(1)]
      g=g(i)
    endfunction

    function [z,lhs]=ystart(x,z,M)
      //z=[4*x*(1-x)*M ; 4*(1-2*x)*M]
      z=[M;0]
      //lhs=[-exp(4*x*(1-x)*M)]
      lhs=0
    endfunction

    for M=[1 4]
       if M==1
          z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,M),tol=tol);
       else
          z1=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,M),tol=tol);
       end
    end

    // Integrating the ode yield e.g. the two solutions yex and yex1.

    function y=f(c),y=c.*(1-tanh(sqrt(c)/4).^2)-2;endfunction
    c=fsolve(2,f);

    function y=yex(x,c)
      y=log(c/2*(1-tanh(sqrt(c)*(1/4-x/2)).^2))
    endfunction

    function y=f1(c1), y=2*c1^2+tanh(1/4/c1)^2-1;endfunction
    c1=fsolve(0.1,f1);

    function y=yex1(x,c1)
      y=log((1-tanh((2*x-1)/4/c1).^2)/2/c1/c1)
    endfunction

    disp(norm(z(1,:)-yex(x)),'norm(yex(x)-z(1,:))= ')
    disp(norm(z1(1,:)-yex1(x)),'norm(yex1(x)-z1(1,:))= ')

    clf();
    subplot(2,1,1)
    plot2d(x,z(1,:),style=[5])
    xtitle('Two different solutions','x',' ')
    subplot(2,1,2)
    plot2d(x,z1(1,:),style=[5])
    xtitle(' ','x',' ')

    // example #4. A multi-point boundary value problem.

    // DE y'''(x)=1
    // z=[y(x);y'(x);y''(x)]
    // BV: y(-1)=2 y(1)=2
    // Side condition: y(0)=1

    a=-1;b=1;c=0;
    // The side condition point c must be included in the array fixpnt.
    n=1;
    m=[3];

    function rhs=fsub(x,z)
      rhs=1
    endfunction

    function g=gsub(i,z)
      g=[z(1)-2 z(1)-1 z(1)-2]
      g=g(i)
    endfunction

    N=10;
    zeta=[a c b];
    x=linspace(a,b,N);

    z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,fixpnt=c);

    function y=yex(x)
    y=x.^3/6+x.^2-x./6+1
    endfunction

    disp(norm(yex(x)-z(1,:)),'norm(yex(x)-z(1,:))= ')

### See Also

`bvode`, `ode`, `dassl`,

### Author

Rainer von Seggern