shooting method for coupled ODE

Question

T K 2020-4-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/519875-shooting-method-for-coupled-ode

回答： T K 2022-8-15

dd
ggg
gg

13 个评论
显示 11更早的评论隐藏 11更早的评论

T K 2020-4-24

在 MATLAB Online 中打开

error 
Undefined function or variable 'dg'.
Error in proj/projfun (line 69)
        dy(2) = (E*f)+(E*x*df)+(f*df)+((dg)/(alfa*Pr*rho))-(df/x)-(f/(x*x));
Error in bvparguments (line 105)
    testODE = ode(x1,y1,odeExtras{:});
Error in bvp4c (line 130)
    bvparguments(solver_name,ode,bc,solinit,options,varargin);
Error in proj (line 37)
    sol= bvp4c(@projfun,@projbc,solinit,options);
 
 (general code is attached as m file )
function sol= proj
clc;clf;clear;
%y1=F, y2=F', y3=t, y4=dt, y5=C, y6=dC, y7=E,y8=E',y9=g
KB        = 1.3807e-15;%Boltzman Constant
Mu        = 1e-2;
K         = 1e5;
Kp        = 40e5;
dp        = (100e-9)*1e2;
rho       = 997.1/1000;
C         = 4179e4;
alfa      = K/(rho*C);
%rhonf     = 3970/1000;
%cnf      = 765e4;
B         = KB/(3*pi*Mu*dp);
Bs        = ((0.26*K)/(2*K+Kp))*(Mu/rho);
myLegend1 = {};
myLegend2 = {};
rr = [4 5 6 7];
%pp =[0.3 0.6 0.7 0.9];
%qq =[4.81 4.9 4.95 5];
for i = 1:numel(rr)
   Pr = rr(i);
%    gamma=pp(i);
    %gammmafi=qq(i);
    %Pr=7;
  Tinf      = 300;  
 Fiinf     = 0.045;
gamma     = 0.05546;
gammafi   = 0.0077;
DT        = gamma*Tinf;
Dfi       = gammafi*Fiinf;
   y0 = [1,0,1,0,1,0,1,0,1];
options =bvpset('stats','on','RelTol',1e-4);
%m = linspace(0.01,1);
    m = linspace(0.1,0.44);
solinit = bvpinit(m,y0);
    sol= bvp4c(@projfun,@projbc,solinit,options);
    %y1=deval(sol,0)
    solinit= sol;
    figure(1)
    plot(sol.x,(sol.y(9,:)))
    grid on,hold on
    myLegend1{i}=['pr = ',num2str(rr(i))];
    figure(2)
    plot(sol.x,sol.y(9,:))
    grid on,hold on
    myLegend2{i}=['Pr= ',num2str(rr(i))];
    i=i+1;
end
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)
    function dy= projfun(x,y)
        dy= zeros(9,1);
%         alignComments
        f  = y(1);
        df = y(2);
        t  = y(3);
        dt = y(4);
        c  = y(5);
        dc = y(6);
        E  = y(7);
        dE = y(8);
        g  = y(9);
        dy(1) = df;
        dy(2) = (E*f)+(E*x*df)+(f*df)+((dg)/(alfa*Pr*rho))-(df/x)-(f/(x*x));
        dy(3) = dt;
        dy(4) = (Pr*dt*x*E)+(Pr*f*dt)-(dt/x);
        dy(5) = dc;
        dy(6) = (1/(B*Tinf))*((alfa*Pr*dc*f)+(alfa*Pr*E*x*dc)-(((DT*Bs*Fiinf)/(Tinf*Dfi))*((dt/x)+((Pr*dt*x*E)+(Pr*f*dt)-(dt/x)))))-(dc/x);
        dy(7) = (-1/x)*(4*E+df+(f/x));
        dy(8) = (4*E*E) +(E*x*((-1/x)*(4*E+df+(f/x)))+f*((-1/x)*(4*E+df+(f/x))) -(((-1/x)*(4*E+df+(f/x)))/x));
        dy(9) = dg ;
    end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
    ya(3)-1; 
    ya(5)-1;
    ya(7)-1;
    ya(9)-1;
    yb(3);
    yb(5);
    yb(7);
    yb(8)];
end

T K 2020-4-26

在 MATLAB Online 中打开

function sol= proj
clc;clf;clear;
%y1=F, y2=F', y3=t, y4=dt, y5=C, y6=dC, y7=E,y8=E',y9=g
KB        = 1.3807e-15;%Boltzman Constant
Mu        = 1e-2;
K         = 1e5;
Kp        = 40e5;
dp        = (100e-9)*1e2;
rho       = 997.1/1000;
C         = 4179e4;
alfa      = K/(rho*C);
%rhonf     = 3970/1000;
%cnf      = 765e4;
B         = KB/(3*pi*Mu*dp);
Bs        = ((0.26*K)/(2*K+Kp))*(Mu/rho);
myLegend1 = {};
myLegend2 = {};
rr = [4 5 6 7];
%pp =[0.3 0.6 0.7 0.9];
%qq =[4.81 4.9 4.95 5];
for i = 1:numel(rr)
   Pr = rr(i);
%    gamma=pp(i);
    %gammmafi=qq(i);
    %Pr=7;
  Tinf      = 300;  
 Fiinf     = 0.045;
gamma     = 0.05546;
gammafi   = 0.0077;
DT        = gamma*Tinf;
Dfi       = gammafi*Fiinf;
   y0 = [1,0,1,0,1,0,1,0,1,0];
options =bvpset('stats','on','RelTol',1e-4);
%m = linspace(0.01,1);
    m = linspace(0.008,0.2);
solinit = bvpinit(m,y0);
    sol= bvp4c(@projfun,@projbc,solinit,options);
    %y1=deval(sol,0)
    solinit= sol;
    figure(1)
    plot(sol.x,(sol.y(9,:)))
    grid on,hold on
    myLegend1{i}=['pr = ',num2str(rr(i))];
    figure(2)
    plot(sol.x,sol.y(10,:))
    grid on,hold on
    myLegend2{i}=['Pr= ',num2str(rr(i))];
    i=i+1;
end
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)
    function dy= projfun(x,y)
        dy= zeros(9,1);
%         alignComments
        f  = y(1);
        df = y(2);
        t  = y(3);
        dt = y(4);
        c  = y(5);
        dc = y(6);
        E  = y(7);
        dE = y(8);
        g  = y(9);
        dg=y(10);
        dy(1) = df;
        dy(2) = (E*f)+(E*x*df)+(f*df)+((dg)/(alfa*Pr*rho))-(df/x)-(f/(x*x));
        dy(3) = dt;
        dy(4) = (Pr*dt*x*E)+(Pr*f*dt)-(dt/x);
        dy(5) = dc;
        dy(6) = (1/(B*Tinf))*((alfa*Pr*dc*f)+(alfa*Pr*E*x*dc)-(((DT*Bs*Fiinf)/(Tinf*Dfi))*((dt/x)+((Pr*dt*x*E)+(Pr*f*dt)-(dt/x)))))-(dc/x);
        dy(7) = (-1/x)*(4*E+df+(f/x));
        dy(8) = (4*E*E) +(E*x*((-1/x)*(4*E+df+(f/x)))+f*((-1/x)*(4*E+df+(f/x))) -(((-1/x)*(4*E+df+(f/x)))/x));
        dy(9) = dg ;
        dy(10)=0;
    end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
    ya(3)-1; 
    ya(5)-1;
    ya(7)-1;
    ya(9)-1;
    yb(1)-0.1;
    yb(3);
    yb(5);
    yb(7);
    yb(8)];
end

T K 2020-4-27

https://www.mathworks.com/matlabcentral/answers/519875-shooting-method-for-coupled-ode#comment_835603

J. Alex Lee 2020-4-27

在 MATLAB Online 中打开

You've asserted that g' does not vary in x by the line, i.e., g''=0

dy(10)=0;

which may or may not be true of your equations, I don't know.

So rather than g being constant (as I suggested before incorrectly), g' will be constant and take whatever value is consistent with your new boundary condition

yb(1)-0.1

Since I don't know the origin of this BC, I don't know if it is consistent with your original problem. I'm not used to dealing with implicit forms, maybe I'm missing something basic. Hope someone else can chime in at this point, I'm at my limit.

请先登录，再进行评论。

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