shooting method for coupled ODE

dd
ggg
gg

[J. Alex Lee]

In your code, you're setting E' to g' in line 83...but really you want the RHS of line 84 to be assigned to dy(7)...right?

Is your issue that you have a BC on g itself, but you aren't solving for g?

[J. Alex Lee]

Please wrap your updated code in the code wrapper (use alt+enter).

What error are you seeing?

[T K]

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]

Doctor J. Alex Alee, how are you? Is it new for the sent code ??

[J. Alex Lee]

The error is straightforward...dg is undefined.

The reason is that you don't have equations in a form where all of the derivatives are explicitly written out (you don't have an explicit form dg/dx = ...)

It's not clear to me that it is possible to express your system in terms of 9 explicit derivatives.

If not, here's somet things to get started:

https://www.mathworks.com/matlabcentral/answers/99336-can-an-ode-be-solved-with-an-implicit-form-in-matlab

https://www.mathworks.com/help/matlab/ref/ode15i.html

I don't think you would be able to use the "shooting" method either.

[T K]

Welcome Doctor J.Alex Lee .

I am very happy and grateful for your effort with me.

I have expressed my system with 10 explicit derivatives As shown in line (77&78) in the code and there is no error appear when run the code.

From your point of view, is the code correct (by expressing the system with 10 explicit derivatives instead of 9 explicit derivatives) ????

[T K]

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

[J. Alex Lee]

You are now just asserting that , i.e., that g is constant (with value 1 according to your 5th BC). It's unlikely what you want. If it is, just remove g completely from the equations and insert 0 in place of dg in Eq. 2.

I'm also unclear where the new BC is coming from (yb(1)-0.1=0) associated with the new DoF.

So no, I don't think this is what you want...

[T K]

Doctor J. Alex, sorry for Bc of yb(1), will ignore it. The question, Doctor, if I ignore g of equations, how do I draw them ?? I want to draw it (g is the pressure)

[J. Alex Lee]

I do not understand what you mean.

[T K]

Excuse me doctor can you help me What are the notes that you mentioned yesterday in this code, especially line 77 & 78?

[T K]

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

[J. Alex Lee]

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.

Answer 1

[T K]

ff

ff