2016-05-17 8 views
0

区分的に定義された関数でbvp4cを使用する際に問題があります。 私はコードをテストし、区分的に定義された関数が一定の場合はうまく動作します。 問題は、区分的に定義された関数が定数ではない領域でグラフに間違った結果が得られることです(確かに分かります)。Matlabのbvp4cソルバに区分的に定義された関数を含める方法

この問題の対処方法に関するご意見やご提案はありますか?

おかげ

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 

答えて

0

は、汝は突然のジャンプを使用して脆弱なソルバーを驚かせることなかれ。

任意の次数pソルバは、少なくともp回のODE関数が連続的に微分可能であり、グリッドポイントのメッシュを適切に適応させることを期待しています。ローカルステップサイズ偏差があれば、特異点の近くで過度に振動する可能性があり、計算時間が長くなり、ステップサイズがアンダーフローします。

この問題を回避するには、イベント(BVPがサポートされている場合)を使用してモデル/ ODE関数を切り替えるか、提供されたmultipoint mechanismを使用して、次に、多くの分岐がある関数ではなく、単純な配列をパラメータとして使用することもできます。

関連する問題