coupled ode for 2nd order
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaire
Torsten
le 7 Fév 2023
Please put your equations, initial and boundary conditions into a readable format.
Réponses (1)
Torsten
le 7 Fév 2023
Modifié(e) : Torsten
le 7 Fév 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 commentaires
Torsten
le 7 Fév 2023
Modifié(e) : Torsten
le 7 Fév 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.
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!