A continuous angular function defined over 3 intervals [0, theta2], [theta2, theta1] and [theta1, pi] can't compile/ show continuity (although it is theoretically).

Question

Mahmoud Zguidi 2023-5-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1956309-a-continuous-angular-function-defined-over-3-intervals-0-theta2-theta2-theta1-and-theta1-pi

编辑： Dyuman Joshi 2023-9-12

在 MATLAB Online 中打开

The function is briefly described as in the attached images:

The following are the so-called boundary conditions:

The following is a continuity relation between theta1 andc theta2:

I would like to plot the function theta(z).

These are 2 slightly different codes for the same problem:

%CODE1:
clear;
clc;
% Define parameters
delta_b1 = 0.05;
delta_b2 = 0.1;
h_1 = 0.2;
h_2 = 0.05;
d = 0.01;
theta_2 = pi/3;
% Define integrands
A = @(theta_i) delta_b1./(d.*sqrt(abs((sin(theta_i)).^2 + 2.*h_1.*(1-cos(theta_i)))));
B = @(theta_i) delta_b1./(d.*sqrt(abs((sin(theta_i)).^2 - 2.*h_1.*(1+cos(theta_i)))));
eqn = @(theta_i) -(2.*h_2./delta_b2.^2).*cos(theta_i) + (delta_b2.^(-2)).*(sin(theta_i)).^2 + ...
    ((delta_b1.^(-2))-(delta_b2.^(-2))).*((sin(theta_2)).^2 - (sin(theta_i)).^2)...
    - 2.*(h_1./delta_b1.^2 - h_2./delta_b2.^2).*(cos(theta_2) - cos(theta_i))...
    + 2.*h_1./delta_b1^2;
C = @(theta_i) 1./(d.*sqrt(abs(eqn(theta_i) + 2.*h_1./delta_b1^2)));
% Find theta1 and theta2
syms x
eqn_sym = eqn(x);
theta_1 = vpasolve(eqn_sym, x, [0, pi]);
theta_1 = double(theta_1);
theta_2 = vpasolve(eqn_sym, x, [theta_1, pi]);
theta_2 = double(theta_2);
% Integrate the integrands
theta1_int = integral(A, 0, theta_1);
theta2_int = integral(B, theta_1, theta_2);
theta3_int = integral(C, theta_2, pi);
% Compute theta_g
theta_g = (theta1_int + theta2_int + theta3_int)/pi;
% Plot the function
theta_range = linspace(0, pi, 200);
theta_vals = arrayfun(@(theta_i) piecewise(theta_i<theta_1, A(theta_i), theta_i>theta_2, B(theta_i), C(theta_i)), theta_range);
plot(theta_range, theta_vals, 'LineWidth', 2);
% Label and format the plot
xlabel('$\xi$', 'interpreter', 'latex')
ylabel('$\theta(\xi)$', 'interpreter', 'latex')
title({'Numerical simulation of the polar angle $\theta(\xi)$ along $\xi$ interval $[0, \pi]$'}, 'interpreter', 'latex');
legend('interpreter', 'latex', 'Location', 'northeast', 'NumColumns', 1, 'FontSize', 12);
grid on;
ax = gca;
ax.TickLabelInterpreter = 'latex';
%%%
%%%
%CODE2:
clear;
clc;
figure
grid on
hold on
theta_2 = pi/3;
delta_b1 = 3;
delta_b2 = 2;
d = 5;
h_1 = 0.25;
h_2 = 0.125;
syms x
eqn = ((delta_b1.^(-2))-(delta_b2.^(-2))).* ((sin(theta_2)).^2 - (sin(x)).^2)...
    - 2.*(h_1./delta_b1.^2 - h_2./delta_b2.^2).* (cos(theta_2) - cos(x))...
    + 2.*h_1./delta_b1^2 == 0;
S = solve(eqn, x);
theta_1 = abs(double(S(1)))
% The integrand for interval [d, +inf[
A = @(theta_i , delta_b1, h_1, d) delta_b1./(d.* sqrt(abs((sin(theta_i )).^2 + 2.*h_1.*(1-cos(theta_i)))));
% The integrand for interval ]-inf -d]
B = @(theta_i , delta_b1, h_1, d) delta_b1./ (d.* sqrt(abs((sin(theta_i )).^2 -2.*h_1.*(1+cos(theta_i)))));
% The integrand for interval [-d,d]
C = @(theta_i , delta_b1, h_1, delta_b2, h_2, d) 1./(d.* sqrt(abs(-(2.* h_2./ delta_b2.^2) .* cos(theta_i)...
    + (delta_b2.^(-2)).* (sin(theta_i)).^2 + ((delta_b1.^(-2))-(delta_b2.^(-2))).* (sin(theta_2)).^2 ...
    - 2.*((h_1./(delta_b1.^(2)))- (h_2./(delta_b2.^(2)))).* cos(theta_2) + 2.* h_1./delta_b1^.2)));
% the z values
z_1 = @(q1)arrayfun(@(theta)quadgk(@(theta_i) A(theta_i , delta_b1, h_1, d), theta, theta_2), q1);
z_2 = @(q1)arrayfun(@(theta)quadgk(@(theta_i) B(theta_i , delta_b1, h_1, d), theta_1, theta), q1);
z_3 = @(q1)arrayfun(@(theta)quadgk(@(theta_i) C(theta_i , delta_b1, h_1, delta_b2, h_2, d), theta, theta_1), q1);
z_4 = @(q1)arrayfun(@(theta)quadgk(@(theta_i) C(theta_i , delta_b1, h_1, delta_b2, h_2, d), theta_2, theta), q1);
%plots
theta = linspace(0, theta_2);
a1 = z_1(theta)  + 1;
b = theta*180/pi;
plot (a1, b,'Displayname', '$ \theta (\xi) for z \in [d, +\infty[$' , 'LineWidth', 2)
theta = linspace(theta_1, pi);
a2 = -z_2(theta) - 1;
b = theta*180/pi;
plot (a2, b,'Displayname', '$ \theta (\xi) for z \in ]-\infty,-d]$' , 'LineWidth', 2)
theta = linspace(theta_2, theta_1);
a3 = z_3(theta) - 1;
a4 = -z_4(theta) + 1;
b = theta*180/pi;
plot (a3, b,'Displayname', '$ \theta (\xi) for z \in [-d, z_0]$' , 'LineWidth', 2)
plot (a4, b, 'Displayname', '$ \theta (\xi) for z \in [z_0, d]$' , 'LineWidth', 2)
legend('interpreter', 'latex', 'Location', 'northeast', 'NumColumns',1, 'FontSize',7)
xlabel('$\xi$', 'interpreter','latex')
ylabel('$\theta(\xi)$', 'interpreter', 'latex')
ax.TickLabelInterpreter = 'latex';
title({' numerical simulation of the polar angle \theta(\xi) along \theta interval [0, \pi]'});
ax = gca;
ax.TitleFontSizeMultiplier = 1;

The plot of this code provides a discontinuous function although the quadgk function is used for maximum numerical precision.

It seems no matter how hard i try the plot isn't good enough. i would appreciate a solution for this small problem. Thank you :)^^

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

SANKALP DEV 2023-9-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1956309-a-continuous-angular-function-defined-over-3-intervals-0-theta2-theta2-theta1-and-theta1-pi#answer_1307866

在 MATLAB Online 中打开

Hello Mahmoud.

I understand that you are currently working on plotting the polar angle 'theta' as a function of 'x.'

After reviewing your original code, it came to my attention that you've utilized numerical integration methods for computing the integrals. However, it's worth noting that this approach may lead to somewhat jagged plots due to the inherent nature of numerical approximations. To achieve smoother plots, I'd like to suggest employing the 'quadgk' function available in MATLAB.

Please consider following code

% Evaluating the functions for the given range
a1 = arrayfun(@(theta_i) quadgk(A, theta_i, theta_2), theta_range) +1 ;
a2 = arrayfun(@(theta_i) -quadgk(B, theta_1, theta_i), theta_range)-1;
a3 = arrayfun(@(theta_i) quadgk(C, theta_i, theta_1), theta_range) -1 ;
a4 = arrayfun(@(theta_i) -quadgk(C, theta_2, theta_i), theta_range) +1;