error while solving the coupled ode

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!

[Torsten]

We cannot tell without testing, and even with testing it's not certain that we can help. So please provide your code.

[KM]

Thanks @Torsten for the comement. and v sorry I forgot to attathed the code.

Here,

[Torsten]

What is "a_tilde" compared to "a" ?

[KM]

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".

Answer 1

[Torsten]

% 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

[KM]

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]

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]

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]

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]

I can not do more than you could do - trial and error.

[KM]

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]

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]

Yes, it's only because of the spherical coordinates.

Basically here, a field is taking shperical ansatz as a solution for field equations.

[Torsten]

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]

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]

These equations look a lot better than the ones you posted first.

Did you try to solve them ?

[KM]

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.