coupled ode for 2nd order

1 次查看（过去 30 天）
KM
KM 2023-2-7
编辑： Torsten 2023-2-7
I have coupled differential equations:
a(r) \sin ^{4} f+2 \cot (f) a_{r} f_{r}+\frac{a_{r}}{r}-a_{r r}=0,{l}\lambda_{0} r(1-\cos f) \sin f\left(\cos f-\cos ^{2} f+\sin ^{2} f\right) \\ \quad+\frac{a^{2} \cos f \sin f}{r}-\frac{\cot f \csc ^{2} f a_{r}^{2}}{r}-f_{r}-r f_{r r}=0
with certain boundary conditions.
  1 个评论
Torsten
Torsten 2023-2-7
Please put your equations, initial and boundary conditions into a readable format.

请先登录，再进行评论。

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

回答（1 个）

Torsten
Torsten 2023-2-7
编辑：Torsten 2023-2-7
Ran in:
I didn't compare equations and boundary conditions with those listed above.
% 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 = 10; % Number of initial grid points used by bvpc4
tol = 1e-3; % Maximum allowed relative error.
L = 10;
N = 2;
n = 0;
m = 0;
g = 5;
lambda = 0;
% Boundary conditions
y0 = [0, -1, N*pi, 0];
% Initial conditions
A = @(xi) (1-tanh(((L*xi)/R)-(L/3)))/2;
dA = cosh(theta)*(coth(delta)-delta*csch(delta)^2);
F = @(xi) (1+tanh(((L*xi)/R)-(L/3)))/2;
dF = (1-delta*coth(delta))*csch(delta);
solinit = bvpinit(linspace(delta, R, initialPoints), [A(delta), F(delta), dA,dF]);
% 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(@(xi, y) heatGauge(xi, y, lambda, g, m, n), @(ya, yb) bcheatGauge(ya, yb, y0),...
solinit, options);
xi = linspace(delta, R, 1e4);
y = deval(sol, xi);
plot(xi,y)
function dy1 = heatGauge(xi, y, lambda, g, m, n)
dy1 = [y(3)...
y(4)...
y(3)./xi + (g^2 * (1+y(1)) * (1+(lambda^2*y(4)^2)) * sin(y(2))^2)...
(1./(1+(lambda^2*(n*(y(1)+1)./xi).^2*sin(y(2))^2))) .* ( ((sin(y(2))*cos(y(2))*(n*(y(1)+1)./xi).^2) + (m^2*sin(y(2)))) - (y(4)./xi).*( ((lambda^2*(n*(y(1)+1)./xi)*sin(y(2))^2).*((n*(y(1)+1)./xi)+(2*xi.*(((xi*n*y(3))-(n*y(1))-n)./xi.^2)))) + 1 + (lambda^2*y(4)*xi.*(n*(y(1)+1)./xi).^2.*sin(y(2))*cos(y(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
  3 个评论
显示 1更早的评论
KM
KM 2023-2-7
lambda_0 is just a constant.
Torsten
Torsten 2023-2-7
编辑：Torsten 2023-2-7
Then you should remember what changes you made to the function because this one worked:
dy1 = [y(3)...
y(4)...
y(3)./xi + (g^2 * (1+y(1)) * (1+(lambda^2*y(4)^2)) * sin(y(2))^2)...
(1./(1+(lambda^2*(n*(y(1)+1)./xi).^2*sin(y(2))^2))) .* ( ((sin(y(2))*cos(y(2))*(n*(y(1)+1)./xi).^2) + (m^2*sin(y(2)))) - (y(4)./xi).*( ((lambda^2*(n*(y(1)+1)./xi)*sin(y(2))^2).*((n*(y(1)+1)./xi)+(2*xi.*(((xi*n*y(3))-(n*y(1))-n)./xi.^2)))) + 1 + (lambda^2*y(4)*xi.*(n*(y(1)+1)./xi).^2.*sin(y(2))*cos(y(2))) ) ) ];
Maybe because the ... are missing in the third line ?
And remember that in your new code, you use cot(x) which is Inf at all multiples of pi. This can easily lead to a singular Jacobian.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by