coupled ode for 2nd order
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentar
Torsten
am 7 Feb. 2023
Please put your equations, initial and boundary conditions into a readable format.
Antworten (1)
Torsten
am 7 Feb. 2023
Bearbeitet: Torsten
am 7 Feb. 2023
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 Kommentare
Torsten
am 7 Feb. 2023
Bearbeitet: Torsten
am 7 Feb. 2023
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.
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!