Need the graph for F,G,theta

conedisk2()

function conedisk2
    % Parameter values
    A1 = 1.10629; 
    A2 = 1.15;
    A3 = 1.15;
    A4 = 1.15;
    A5 = 1.15;
    A6 = 1.15;
    M = 0.2;
    Grt = 5;
    Pr = 0.71;
    R = 0.2;
    Ec = 0.1;
    Q = 0.1;
    Rew = 12;
    Red = 12;
    n = -1;
    g1 = 0;
    g2 = 0;
    g3 = 0;
    g4 = 0;
    inf = 0.06992681;
    
    % Set solver options to increase the maximum number of mesh points
    options = bvpset('RelTol', 1e-5, 'AbsTol', 1e-7, 'NMax', 5000);
    
    % Defining parameters 
    solinit = bvpinit(linspace(0, 1, 200), [0 g1 g2 Rew g3 0 1 g4]);
    sol = bvp4c(@bvp2D, @bc2D, solinit);
    x = sol.x;
    y = sol.y;
    
    % Plotting of the velocity
    figure(1)
    plot(x, y(1,:), 'linewidth', 1)
    hold on
    xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('F(\eta)', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    % Plotting of G
    figure(2)
    plot(x, y(5,:), 'linewidth', 1)
    hold on
    xlabel('\eta ', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('G(\eta) ', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    % Plotting of theta
    figure(3)
    plot(x, y(7,:), 'linewidth', 1)
    hold on
    xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('\theta(\eta)', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    % Residual of the boundary conditions
    function res = bc2D(y0, yinf)
        res = [y0(1);yinf(1); y0(4); yinf(4);y0(5)-Rew; yinf(5)-Red;y0(7)-1; yinf(7)];
    end
    % System of First Order ODEs
    function dydx = bvp2D(eta, y)
        yy1 = (-(1+eta^2)*((10*A1*eta+A2*eta*y(1)-A2*y(4))*y(3))-3*(2*A1+7*A1*(eta^2)+A2*(2*eta^2+1)*y(1)-A2*eta*y(4))*y(2)-2*A2*y(5)*y(6)-3*(A1+A2*y(1))*y(4)+A3*M*y(2)-A4*Grt*y(8))/(A1*(1+eta^2)^2);
        yy2 = eta*y(2);
        yy3 = (-3*A1*eta*y(6)-A2*(eta*y(1)-y(4))*y(6)+A3*M*y(5))/(A1*(eta^2+1));
        yy4 = ((A5+4/3*R)*((eta*((2*n)-1))*y(8)+n^2*y(7))-A3*Pr*M*Ec*(y(1)^2+y(5)^2)-Pr*Q*y(7)+A6*Pr*(y(4)*y(8)+n*y(1)*y(7)-eta*y(1)*y(8)))/(A5+4/3*R)*(1+eta^2);
        dydx = [y(2); y(3); yy1; yy2; y(6); yy3; y(8); yy4];
    end
end 

The code is running but I am unable to get the graph what i have attached for radial velocity and azimuthal velocity

[Torsten]

The error can lie in your parameters, your equations, your boundary conditions ... How should we be able to help in this case ?

Answer 1

[jinjin lee]

you want get F,G,theta？just output the data。

data=conedisk2()

Function

function data=conedisk2()
    % Parameter values
    A1 = 1.10629; 
    A2 = 1.15;
    A3 = 1.15;
    A4 = 1.15;
    A5 = 1.15;
    A6 = 1.15;
    M = 0.2;
    Grt = 5;
    Pr = 0.71;
    R = 0.2;
    Ec = 0.1;
    Q = 0.1;
    Rew = 12;
    Red = 12;
    n = -1;
    g1 = 0;
    g2 = 0;
    g3 = 0;
    g4 = 0;
    inf = 0.06992681;
    
    % Set solver options to increase the maximum number of mesh points
    options = bvpset('RelTol', 1e-5, 'AbsTol', 1e-7, 'NMax', 5000);
    
    % Defining parameters 
    solinit = bvpinit(linspace(0, 1, 200), [0 g1 g2 Rew g3 0 1 g4]);
    sol = bvp4c(@bvp2D, @bc2D, solinit);
    x = sol.x;
    y = sol.y;
    
    % Output data
    data.x=x;
    data.F=y(1,:);
    data.G=y(5,:);
    data.theta=y(7,:);
    
    % Plotting of the velocity
    figure(1)
    plot(x, y(1,:), 'linewidth', 1)
    hold on
    xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('F(\eta)', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    % Plotting of G
    figure(2)
    plot(x, y(5,:), 'linewidth', 1)
    hold on
    xlabel('\eta ', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('G(\eta) ', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    % Plotting of theta
    figure(3)
    plot(x, y(7,:), 'linewidth', 1)
    hold on
    xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('\theta(\eta)', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    % Residual of the boundary conditions
    function res = bc2D(y0, yinf)
        res = [y0(1);yinf(1); y0(4); yinf(4);y0(5)-Rew; yinf(5)-Red;y0(7)-1; yinf(7)];
    end
    % System of First Order ODEs
    function dydx = bvp2D(eta, y)
        yy1 = (-(1+eta^2)*((10*A1*eta+A2*eta*y(1)-A2*y(4))*y(3))-3*(2*A1+7*A1*(eta^2)+A2*(2*eta^2+1)*y(1)-A2*eta*y(4))*y(2)-2*A2*y(5)*y(6)-3*(A1+A2*y(1))*y(4)+A3*M*y(2)-A4*Grt*y(8))/(A1*(1+eta^2)^2);
        yy2 = eta*y(2);
        yy3 = (-3*A1*eta*y(6)-A2*(eta*y(1)-y(4))*y(6)+A3*M*y(5))/(A1*(eta^2+1));
        yy4 = ((A5+4/3*R)*((eta*((2*n)-1))*y(8)+n^2*y(7))-A3*Pr*M*Ec*(y(1)^2+y(5)^2)-Pr*Q*y(7)+A6*Pr*(y(4)*y(8)+n*y(1)*y(7)-eta*y(1)*y(8)))/(A5+4/3*R)*(1+eta^2);
        dydx = [y(2); y(3); yy1; yy2; y(6); yy3; y(8); yy4];
    end
end