error while solving the coupled ode

2 次查看（过去 30 天）
KM
KM 2023-2-21
评论： KM 2023-2-24
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 个评论
显示 2更早的评论
Torsten
Torsten 2023-2-21
What is "a_tilde" compared to "a" ?
KM
KM 2023-2-21
a_tilde = n(a+1)/r
The origian equation was in a, "a_tilde" was a substitute.
I have written the full equation in code after the substituting the value of "a_tilde".

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2023-2-21
Ran in:
% 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 个评论
显示 10更早的评论
Torsten
Torsten 2023-2-24
These equations look a lot better than the ones you posted first.
Did you try to solve them ?
KM
KM 2023-2-24
Hi @Torsten. I have got the conditional plots
  1. Code doesn't run for R>2 and g>1.4
  2. and have to modity the bc for a, it's not running for -1 but -0.9. (in the codes)
delta = 0.0001; % Lower integral bound
R = 1; % Upper integral bound
theta = 0; % ArcTan(q/g)
maxPoints = 1e4; % Maximum numer of grid point used by bvpc4
initialPoints = 10; % Number of initial grid points used by bvpc4
tol = 1e-4;
L = 10; % Maximum allowed relative error.
g = 0.0;
mu = sqrt(0.1);
n = 1;
lambda = 1;
% Boundary conditions
ya(1) = 0;
yb(1) =-0.9;
ya(2) = 1;
yb(2) = 0;
y0 = [0, -0.9, 1, 0];
% Initial conditions
A = @(xi) 3*xi./sinh(3*xi)-1;
F = @(xi) 3*xi./sinh(3*xi);
dA = (1-delta*coth(delta))*csch(delta);
dF = (1-delta*coth(delta))*csch(delta);
solinit = bvpinit(linspace(delta, R, initialPoints), [A(delta), F(delta), dA, dF]);
options = bvpset('RelTol', tol, 'NMax', maxPoints);
Plot(xi.^2/2 , y(1,:), xi.^2/2, y(2,:))
Although I was looking for plots for at most g = 10 and R = 8.
I can't figure out these any further. If any input from your end can make it as expected, thanks in advance.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by