error is not showing but plot is not generating
조회 수: 1 (최근 30일)
이전 댓글 표시
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
답변 (1개)
Walter Roberson
2023년 12월 2일
편집: Walter Roberson
2023년 12월 4일
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
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.
참고 항목
카테고리
Help Center 및 File Exchange에서 Structures에 대해 자세히 알아보기
제품
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!