error while solving the coupled ode
3 次查看(过去 30 天)
显示 更早的评论
Hello,
I have been trying to solve the following ode
but a "singular jacobian error" is recurring and is very sensitive to parameters. I think it's because of the guess function issue but I'm not very sure of it. I tried many guess functions but didn't seem to work. If there is any way out or suggestion to overshoot this error, please do help.
Thanks!
4 个评论
采纳的回答
Torsten
2023-2-21
% Defining parameters
delta = 0.02; % Lower integral bound
R = 5; % Upper integral bound
theta = 0; % ArcTan(q/g)
maxPoints = 1e6; % Maximum numer of grid point used by bvpc4
initialPoints = 100; % Number of initial grid points used by bvpc4
tol = 1e-3; % Maximum allowed relative error.
L = 10;
N = 1;
n = 1;
m = 0;
g = 5;
lambda = 1;
% Boundary conditions
y0 = [0, -1, N*pi, 0];
% Initial conditions
A = @(r) (1-tanh(((L*r)/R)-(L/3)))/2;
dA = cosh(theta)*(coth(delta)-delta*csch(delta).^2);
F = @(r) (1+tanh(((L*r)/R)-(L/3)))/2;
dF = (1-delta*coth(delta))*csch(delta);
solinit = bvpinit(linspace(delta, R, initialPoints), [1 1 1 1]);
% Solves system using bvpc4
options = bvpset('RelTol', tol, 'NMax', maxPoints); % This function sets the allowed
%relative error and maximum number of grid points.
sol = bvp4c(@(r, y) heatGauge(r, y, lambda, g, m, n), @(ya, yb) bcheatGauge(ya, yb, y0),...
solinit, options);
r = linspace(delta, R, 1e4);
y = deval(sol, r);
plot(r,y(1,:),r, y(2,:))
grid on
function dy = heatGauge(r, y, lambda, g, m, n)
a = y(1);
f = y(2);
adot = y(3);
fdot = y(4);
atilde = n*(a+1.0)/r;
dy(1) = adot;
dy(2) = fdot;
dy(3) = a/r + g^2*(1+a)*(1+lambda^2*fdot^2)*sin(f)^2;
dy(4) = (-fdot/r*((2*n*adot-atilde)*lambda^2*atilde*sin(f)^2+lambda^2*r*fdot*atilde^2*sin(f)*cos(f)+1)...
+atilde^2*sin(f)*cos(f)+m^2*sin(f))/(1+lambda^2*atilde^2*sin(f)^2);
end
function res = bcheatGauge(ya, yb, y0)
res = [ya(1) - y0(1);yb(1) - y0(2);ya(2) - y0(3);yb(2) - y0(4)];
end
12 个评论
Torsten
2023-2-24
These equations look a lot better than the ones you posted first.
Did you try to solve them ?
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!