can it be possible to use odeVectorField

Question

MINATI 2020-5-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/530143-can-it-be-possible-to-use-odevectorfield

评论： MINATI 2020-5-25

syms f(x) g(x) h(x) Pr Nb Nt Le
[V] = odeToVectorField( diff(f,3) - diff(f)^2 + f * diff(f,2) == 0,diff(g,2) + Pr * ( f * diff(g) + Nb * diff(g)* diff(h)+ Nt * diff(g)^2) == 0,...
                                       diff(h,2) + Le * f * diff(h) + ( Nt/Nb ) * Pr * ( f * diff(g) + Nb * diff(g)* diff(h)+ Nt * diff(g)^2)== 0)
%%  I want to incorporate the expression of [V]  into below code, but I am getting wrong expression (not as  dydx  below) 
%%%  can it be used in   bvp 4c  code  and  in which  style
function main
Pr=10; Le=10; Nt=.1;Nb=.1;
xa=0;xb=6;solinit=bvpinit(linspace(xa,xb,1000),[0 1 1 1 0 0 0]);options=[];
sol=bvp4c(@ode,@bc,solinit,options);xint=linspace(xa,xb,100);S=deval(sol,xint);
function res=bc(ya,yb)
res=[ya(1); ya(2)-1; ya(4)-1; ya(6)-1; yb(2); yb(4); yb(6)];        
end 
function dydx=ode(x,y)    
dydx=[y(2);  y(3);  y(2)^2 - y(3)*y(1);    y(5);          - Pr*(y(1)*y(5) + Nb*y(5)*y(7) + Nt*y(5)^2);
           y(7);    - Le*y(7)*y(1) + (Nt/Nb)*Pr*(y(1)*y(5) + Nb*y(5)*y(7) + Nt*y(5)^2)];             
end
end

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

MINATI PATRA 2020-5-23

Dear Darova that M gives wrong code as you can compare and another question is how to include in bvp 4c code

MINATI 2020-5-24

Please anybody invest some time here to help me.

请先登录，再进行评论。

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

Answer 1

darova 2020-5-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/530143-can-it-be-possible-to-use-odevectorfield#answer_437613

Try this

clc,clear
syms f(x) g(x) h(x) x Y
Pr=10;
Le=10;
Nt=.1;
Nb=.1;
[V,y1] = odeToVectorField( diff(f,3) - diff(f)^2 + f * diff(f,2) == 0,...
    diff(g,2) + Pr * ( f * diff(g) + Nb * diff(g)* diff(h)+ Nt * diff(g)^2) == 0,...
    diff(h,2) + Le * f * diff(h) + ( Nt/Nb ) * Pr * ( f * diff(g) + Nb * diff(g)* diff(h)+ Nt * diff(g)^2)== 0)
F = matlabFunction(V,'vars', [x Y]);
bc = @(ya,yb) [ya(1)    % g
    ya(2)-1 % Dg
    ya(4)-1 % f
    ya(6)-1 % Df
    yb(2)   % D2f
    yb(4)   % h
    yb(6)]; % dh
xa=0;
xb=6;
solinit = bvpinit(linspace(xa,xb,1000),[0 1 1 1 0 0 0]);
options = [];
sol = bvp4c(F,bc,solinit,options);
plot(sol.x,sol.y)