error is not showing but plot is not generating
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Respuestas (1)
Walter Roberson
el 2 de Dic. de 2023
Editada: Walter Roberson
el 4 de Dic. de 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 comentarios
Walter Roberson
el 4 de Dic. de 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.
Ver también
Categorías
Más información sobre Structures en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!