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

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 :)^^

Answer 1

[SANKALP DEV]

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;

To know more about ‘quadgk’ MATLAB function, please refer to the following documentation-https://www.mathworks.com/help/releases/R2020b/matlab/ref/quadgk.html?searchHighlight=quadgk&s_tid=doc_srchtitle

Hope this helps