How to include piecewise defined functions in bvp4c solver?

Question

0 개 추천

I am having trouble using the bvp4c with piecewise defined functions. I tested the code and it works fine when the piecewise defined functions are constant. The problem is that I get wrong results in the graph (that I know for sure) in the area where the piecewise defined functions are not constant, but linear equations.

Any ideas or suggestions on how to handle this problem?

Thanks.

function bvp4
xlow=0;
xhigh=0.30;
solinit=bvpinit(linspace(xlow,xhigh,1000),[0 0]);
sol = bvp4c(@bvp4ode,@bvp4bc,solinit);
xint=[xlow:0.0001:xhigh];
Sxint=deval(sol,xint);
Sxint1=abs(sqrt(Sxint));
xint=[xlow:0.0001:xhigh];
plot(xint,Sxint1(1,:),'r')
function dydx = bvp4ode(x,y)
So=0.00125;
s=1.5;
dydx = [y(2);    
    ((G(x)+125*f(x)*y(1)*(1+1/s^2)^0.5-1000*9.81*So*H(x))/(1000*0.5*l(x)*(f(x)/8)^0.5)-y(2)*2*(-2/3*x+0.071+2/3*0.08)*(-2/3)*b(x))/H(x)/H(x)];
function res = bvp4bc(ya,yb)
res = [ya(1);        yb(1)];
function fval = f(x)
if     (x >= 0) && (x <= 0.08)
    fval = 0.0187;
elseif  (x > 0.08) && (x <= 0.17)
    fval = 0.0298;
elseif  (x > 0.17) && (x <= 0.3)
    fval= 0.0408;
end
function Gval = G(xint)
if     (xint >= 0) && (xint <= 0.08)
    Gval = 0.1306;
elseif  (xint > 0.08) && (xint <= 0.17)
    Gval = 0.1306;
elseif  (xint > 0.17) && (xint <= 0.3)
    Gval = -0.0337;
end
function Hval = H(xint)
if     (xint >= 0) && (xint < 0.08)
    Hval = 0.071;
elseif  (xint >= 0.08) && (xint <= 0.17)
    Hval = -2/3*xint+(0.071+2/3*0.08);
elseif  (xint >0.17) && (xint <= 0.3)
     Hval = 0.011;
end
function bval = b(xint)
if     (xint >= 0) && (xint < 0.08)
    bval = 0;
elseif  (xint >= 0.08) && (xint <= 0.17)
    bval = 1;
elseif  (xint > 0.17) && (xint <= 0.3)
    bval= 0;
end
function lval = l(xint)
if     (xint >= 0) && (xint <= 0.08)
    lval = 0.067;
elseif  (xint > 0.08) && (xint <= 0.17)
    lval = 0.134;
elseif  (xint > 0.17) && (xint <= 0.3)
    lval= 1.165;
end