error while solving the coupled ode

1 次查看（过去 30 天）

显示更早的评论

KM 2023-2-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1915910-error-while-solving-the-coupled-ode

评论： KM 2023-2-24

采纳的回答： Torsten

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 2023-2-21

What is "a_tilde" compared to "a" ?

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 2023-2-21

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1915910-error-while-solving-the-coupled-ode#answer_1176485

在 MATLAB Online 中打开

% 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更早的评论

Torsten 2023-2-24

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

Did you try to solve them ?

KM 2023-2-24

在 MATLAB Online 中打开

Hi @Torsten. I have got the conditional plots

Code doesn't run for R>2 and g>1.4
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,:))