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
i defined all the parameters but plot is not generating any mistakes in the code
回答（1 个）
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 个评论
KARUNA BOPPENA
2023-12-4
Walter Roberson
2023-12-4
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
2023-12-4
Walter Roberson
2023-12-4
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Communications Toolbox 的更多信息
产品
2023-12-2
2023-12-4
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
