error is not showing but plot is not generating
1 view (last 30 days)
Show older comments
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
i defined all the parameters but plot is not generating any mistakes in the code
0 Comments
Answers (1)
Walter Roberson
on 2 Dec 2023
Edited: Walter Roberson
on 4 Dec 2023
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
4 Comments
Walter Roberson
on 4 Dec 2023
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.
See Also
Categories
Find more on Structures in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!