Solve ODE for variable domain

Question

EldaEbrithil 2020-9-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/588607-solve-ode-for-variable-domain

评论： Alan Stevens 2020-9-3

采纳的回答： Alan Stevens

在 MATLAB Online 中打开

Hi all

is it possible to solve the same ODE system for 3 adiacent different domain? Do i need something like that?

if x>x1&&x<=x2
  .
  .
  .
elseif x>=x2&&x<x3
    .
    .
    .
elseif x>x3
    .
    .
    .
end

how continuity will be preserved?

Thank you for the help

Regards

5 个评论
显示 3更早的评论隐藏 3更早的评论

EldaEbrithil 2020-9-3

在 MATLAB Online 中打开

function dYnozzledx=nozzlesinglebobb(x,Ynozzle,xstartconica,xfineconica,Ralfa_end,Rbeta_end,alfa_end_rad,beta_end_rad,xc_end,yc_end,Lnozzle_end,raggio_end,...
    f1,f2,f3)
gamma=1.667;
Pt_nozzle=Ynozzle(1);
M_nozzle=Ynozzle(2);
if (x>=0)&&(x<xstartconica)
Anozzle= pi*((Rbeta_end*cos(beta_end_rad)-x)*tan(beta_end_rad))^2;
Perimetro=2*pi*((Rbeta_end*cos(beta_end_rad)-x)*tan(beta_end_rad));
%p=(Pt_nozzle/(1+((gamma-1)/2)*M_nozzle^2)^(gamma/(gamma-1)))
dA_nozzledx=(2*pi*tan(beta_end_rad)^2)*(x-Rbeta_end*cos(beta_end_rad));
dPt_nozzledx=-Pt_nozzle*(Perimetro*f1*gamma*M_nozzle^2)/(2*Anozzle*(Pt_nozzle/(1+((gamma-1)/2)*M_nozzle^2)^(gamma/(gamma-1)))*Anozzle);
dM_nozzledx=M_nozzle*((-(1+((gamma-1)/2)*M_nozzle^2)/(1-M_nozzle^2))*(dA_nozzledx/Anozzle)+...
    ((1+((gamma-1)/2)*M_nozzle^2)/(1-M_nozzle^2))*(gamma*(M_nozzle^2)*f1*Perimetro/(2*Anozzle*(Pt_nozzle/(1+((gamma-1)/2)*M_nozzle^2)^(gamma/(gamma-1)))*Anozzle)));
elseif (x>=xstartconica)&&(x<=xfineconica)
Anozzle= pi*(yc_end-sqrt(raggio_end^2-(x-xc_end)^2))^2;
Perimetro= 2*pi*(yc_end-sqrt(raggio_end^2-(x-xc_end)^2));
dA_nozzledx=(2*pi*(x-xc_end)/(sqrt(raggio_end^2-(x-xc_end)^2)))*(yc_end-sqrt(raggio_end^2-(x-xc_end)^2));
dPt_nozzledx=-Pt_nozzle*(Perimetro*f2*gamma*M_nozzle^2)/(2*Anozzle*(Pt_nozzle/(1+((gamma-1)/2)*M_nozzle^2)^(gamma/(gamma-1)))*Anozzle);
dM_nozzledx=M_nozzle*((-(1+((gamma-1)/2)*M_nozzle^2)/(1-M_nozzle^2))*(dA_nozzledx/Anozzle)+...
    ((1+((gamma-1)/2)*M_nozzle^2)/(1-M_nozzle^2))*(gamma*(M_nozzle^2)*f2*Perimetro/(2*Anozzle*(Pt_nozzle/(1+((gamma-1)/2)*M_nozzle^2)^(gamma/(gamma-1)))*Anozzle)));
elseif x>xfineconica
    
Anozzle= pi*((x-(Lnozzle_end-(Ralfa_end*cos(alfa_end_rad))))*tan(alfa_end_rad))^2;
Perimetro= 2*pi*((x-(Lnozzle_end-(Ralfa_end*cos(alfa_end_rad))))*tan(alfa_end_rad));
dA_nozzledx=(2*pi*tan(alfa_end_rad)^2)*(x-Lnozzle_end+Ralfa_end*cos(alfa_end_rad));
dPt_nozzledx=-Pt_nozzle*(Perimetro*f3*gamma*M_nozzle^2)/(2*Anozzle*(Pt_nozzle/(1+((gamma-1)/2)*M_nozzle^2)^(gamma/(gamma-1)))*Anozzle);
dM_nozzledx=M_nozzle*((-(1+((gamma-1)/2)*M_nozzle^2)/(1-M_nozzle^2))*(dA_nozzledx/Anozzle)+...
    ((1+((gamma-1)/2)*M_nozzle^2)/(1-M_nozzle^2))*(gamma*(M_nozzle^2)*f3*Perimetro/(2*Anozzle*(Pt_nozzle/(1+((gamma-1)/2)*M_nozzle^2)^(gamma/(gamma-1)))*Anozzle)));
end
dYnozzledx=[dPt_nozzledx;dM_nozzledx];
end

i need to subdived the domain because there are 3 different function that described the Section (Anozzle) trend

EldaEbrithil 2020-9-3

as you can see the central part is a conical section

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2020-9-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/588607-solve-ode-for-variable-domain#answer_489334

在 MATLAB Online 中打开

So you would then call your ode with something like:

[x,Y] = ode45(@nozzlesinglebobb, xspan, Y0,[],xstartconica,xfineconica,Ralfa_end,Rbeta_end,alfa_end_rad,beta_end_rad,xc_end,yc_end,Lnozzle_end,raggio_end,f1,f2,f3);

where xspan = [0 x_final] and Y0 = [Pt0; M0] or similar.

Then

Pt = Y(:,1);
M = Y(:,2);

2 个评论
显示无隐藏无

EldaEbrithil 2020-9-3

编辑：EldaEbrithil 2020-9-3

在 MATLAB Online 中打开

I had written something that

[xSolnozzle,YSolnozzle]=ode23(@(x,Y)nozzlesinglebobb(x,Ynozzle,xstartconica,...
    xfineconica,Ralfa_end,Rbeta_end,alfa_end_rad,beta_end_rad,xc_end,yc_end,Lnozzle_end,raggio_end,...
    f1,f2,f3),xRangenozzle,Ynozzle);

is it correct?