How to include piecewise defined functions in bvp4c solver?

5 views (last 30 days)
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

Accepted Answer

Torsten
Torsten on 18 May 2016
Solve your problem as a multipoint boundary value problem:
Select the intermediate points as 0.08 and 0.17.
Best wishes
Torsten.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!