error is not showing but plot is not generating

clear
%% %% current density equation J=dv/dx=0 at x=0 boundary condition
%% u=y(1)
%% v=y(2)
%% du/dx=dy(1)/dx=y(3)
%%d^2u/dx^2=dy(3)/dx=gamma*y(1)/(1+alpha*y(1)
%% dv/dx=dy(2)/dx=y(4)
%% d^2v/dx^2=dy(4)/dx=a*delta_v*y(4)-(2/epsilon)*gamma*y(1)/(1+alpha*y(1)
alpha = 0.1;
gamma = [1, 50, 100, 500, 1000];
epsilon = 1;
a = 1000;

fcn = @(x, y) [y(3); y(4); (gamma * y(1)) / (1 + alpha * y(1)); a * v * y(4) - (2 / epsilon) * (gamma * y(1)) / (1 + alpha * y(1))];
 bc = @(ya, yb) [ya(1)-1; ya(2); yb(3); 0]; %% J=dv/dx=0 at x=0 so yb(4) is 0
guess = @(x) [1; 0; 0; 0]; %% at x=0, u=1, v=0, du/dx=0, dv/dx=0
xmesh = linspace(0, 1, 100);
solinit = bvpinit(xmesh, guess);
for i = 1:numel(gamma)
    sol = bvp4c(fcn, bc, solinit);
    y_plot = deval(sol, xmesh);
    v = y_plot(1, :);
    dv_dx = y_plot(2, :);
    J = dv_dx;
    delta_v_star = v;

    % Plot delta_v_star versus J
    figure;
    plot(delta_v_star, J);
    xlabel(' delta_v_star');
    ylabel('J');
   title('Plot of delta_v_star verses current density (J)');
   legend('\gamma = 1', '\gamma = 50', '\gamma = 100', '\gamma = 500', '\gamma = 1000');
    xlim([1e-4, 1e2]);
    ylim([0, 70]);
end

Unrecognized function or variable 'v'.

Error in solution>@(x,y)[y(3);y(4);(gamma*y(1))/(1+alpha*y(1));a*v*y(4)-(2/epsilon)*(gamma*y(1))/(1+alpha*y(1))] (line 14)
fcn = @(x, y) [y(3); y(4); (gamma * y(1)) / (1 + alpha * y(1)); a * v * y(4) - (2 / epsilon) * (gamma * y(1)) / (1 + alpha * y(1))];

Error in bvparguments (line 96)
testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 119)
bvparguments(solver_name,ode,bc,solinit,options,varargin);

i defined all the parameters but plot is not generating any mistakes in the code

Answer 1

[Walter Roberson]

The Jacobian is singular if y(4) = dv/dx = 0, which it is because of your boundary conditions.

That said, even if I change the boundary conditions, I still get told singular jacobian.

clear
%% %% current density equation J=dv/dx=0 at x=0 boundary condition
%% u=y(1)
%% v=y(2)
%% du/dx=dy(1)/dx=y(3)
%%d^2u/dx^2=dy(3)/dx=gamma*y(1)/(1+alpha*y(1)
%% dv/dx=dy(2)/dx=y(4)
%% d^2v/dx^2=dy(4)/dx=a*delta_v*y(4)-(2/epsilon)*gamma*y(1)/(1+alpha*y(1)
alpha = 0.1;
gamma = [1, 50, 100, 500, 1000];
epsilon = 1;
a = 1000;

 bc = @(ya, yb) [ya(1)-1; ya(2); yb(3); 0]; %% J=dv/dx=0 at x=0 so yb(4) is 0
guess = @(x) [1; 0; 0; 0]; %% at x=0, u=1, v=0, du/dx=0, dv/dx=0
xmesh = linspace(0, 1, 100);
solinit = bvpinit(xmesh, guess);
syms X Y [1 4]
for i = 1:numel(gamma)
    fcn = @(x, y) [y(3); y(4); (gamma(i) * y(1)) / (1 + alpha * y(1)); a * y(2) * y(4) - (2 / epsilon) * (gamma(i) * y(1)) / (1 + alpha * y(1))];
    i
    F = fcn(X, Y)
    Jac = jacobian(F)
    rank(Jac)
    detJac = det(Jac)
    [N,D] = numden(detJac)
    solve(N)
    solve(D)
    
    sol = bvp4c(fcn, bc, solinit);
    y_plot = deval(sol, xmesh);
    v = y_plot(1, :);
    dv_dx = y_plot(2, :);
    J = dv_dx;
    delta_v_star = v;

    % Plot delta_v_star versus J
    figure;
    plot(delta_v_star, J);
    xlabel(' delta_v_star');
    ylabel('J');
   title('Plot of delta_v_star verses current density (J)');
   legend('\gamma = 1', '\gamma = 50', '\gamma = 100', '\gamma = 500', '\gamma = 1000');
    xlim([1e-4, 1e2]);
    ylim([0, 70]);
end

i = 1

F = 

Jac = 

ans = 4

detJac = 

N = 

D = 

ans = 

0

ans = 

Error using bvp4c
Unable to solve the collocation equations -- a singular Jacobian encountered.

[KARUNA BOPPENA]

Thank you for your response. how could i solve this problem any suggestion please.

[Walter Roberson]

At the very least your y(4) must not start at 0.

However, when I tested with other values, I still got that error, and I am not sure why so.

[KARUNA BOPPENA]

Thank you

[Walter Roberson]

 bc = @(ya, yb) [ya(1)-1; ya(2); yb(3); 0]; %% J=dv/dx=0 at x=0 so yb(4) is 0

That has a fixed element of 0, and expects to work with vectors of length 4. The result can be at most rank 3, so the jacobian of those boundary conditions must be singular.