error while solving the coupled ode

KM
2023 2 21
1 个回答
回答已采纳
更新时间：2023 2 24
2 次查看（30 天）

0 个投票

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更早的评论 隐藏 2更早的评论

Torsten
Torsten 2023-2-21
We cannot tell without testing, and even with testing it's not certain that we can help. So please provide your code.
KM
KM 2023-2-21
Thanks @Torsten for the comement. and v sorry I forgot to attathed the code.
Here,
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:

1 个投票

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更早的评论 隐藏 10更早的评论

KM
KM 2023-2-21
Hi, @Torsten much thanks!
Now the variation like if I increase the range of R to 20, or increase the value of g to 20, jacob error pops up in each case. This error is very sesitive to such variations. How should i play with such error. Also, why the guess function is [1,1,1,1]?
Torsten
Torsten 2023-2-21
Now the variation like if I increase the range of R to 20, or increase the value of g to 20, jacob error pops up in each case. This error is very sesitive to such variations. How should i play with such error.
Trial-And-Error as I did with Tolerances, number of grid points, initial guesses,...
And checking equations for correctness and results for plausibility.
Also, why the guess function is [1,1,1,1]?
Your guess didn't work. So why not [1,1,1,1] ?
KM
KM 2023-2-21
Okay. Thanks!
Your guess didn't work. So why not [1,1,1,1] ?
The guess function is defined considering the boundary conditions of eqns. I have taken the above guess function because they satisfied the nature of curve I was expecting provided they match the boundary conditions as well.
If guess function in a constant function, the variation won't taken place on the cells for the exactness of the solution I expect so that's why.
Torsten
Torsten 2023-2-21
编辑：Torsten 2023-2-21
A guess is only a guess - it won't influence the final solution if the solver is able to converge. It can only help to converge faster or to converge at all.
The guess function is defined considering the boundary conditions of eqns. I have taken the above guess function because they satisfied the nature of curve I was expecting provided they match the boundary conditions as well.
You also supply constant values for all four variables as initial guesses.
[A(delta) F(delta) dA dF]
are all scalars. Function curves are not supplied.
Torsten
Torsten 2023-2-23
I can not do more than you could do - trial and error.
KM
KM 2023-2-23
编辑：KM 2023-2-23
yes, I understand. Thank you!
This particular equation is not giving any solution at all. I tried all parameters I could but always the jacobian error or a mess error.
Torsten
Torsten 2023-2-23
The division by r resp. r^2 : where does it come from ?
Is it because you work in a cylindrical or a spherical coordinate system ?
KM
KM 2023-2-23
Yes, it's only because of the spherical coordinates.
Basically here, a field is taking shperical ansatz as a solution for field equations.
Torsten
Torsten 2023-2-23
Usually, in cylindrical or spherical coordinates, it's not possible to prescribe function values at r = 0.
Thus I doubt that ya(1) and ya(2) can be prescribed. Instead, ya(3) = 0 and ya(4) = 0 are common here.
KM
KM 2023-2-24
编辑：KM 2023-2-24
it's not possible to prescribe function values at r = 0.
Indeed @Torsten and that's why 'r' (=delta) here does not start from 0 but 0.0 in the code to avoid any singularity, but even if I increase delta to 0.1, it still does give singulairy (if that's what jacobian singularity here is).
Thus I doubt that ya(1) and ya(2) can be prescribed.
These are the boundary conditons in terms of spherical coordinates. What I can try it to exclude the computation when part.
Alternatively, I tried below another proposed method but this ain't working to me either.
I can redefine these f and a fields for computations like:
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 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Noncentral t Distribution 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by