Month: December 2025
How to configure the default message and the hyperlink of the help command for published functions
I have published the documentation for a function that I coded with publishing markups in html format but when I type help myFunction in my Command Window at the end appears two lines by default like the following ones:
"Published output in the Help browser
showdemo myFunction (<- hyperlink)."
Is there a way to configure the way this two lines are shown? Like for example, change them to something like this:
"For more information see myFunction"I have published the documentation for a function that I coded with publishing markups in html format but when I type help myFunction in my Command Window at the end appears two lines by default like the following ones:
"Published output in the Help browser
showdemo myFunction (<- hyperlink)."
Is there a way to configure the way this two lines are shown? Like for example, change them to something like this:
"For more information see myFunction" I have published the documentation for a function that I coded with publishing markups in html format but when I type help myFunction in my Command Window at the end appears two lines by default like the following ones:
"Published output in the Help browser
showdemo myFunction (<- hyperlink)."
Is there a way to configure the way this two lines are shown? Like for example, change them to something like this:
"For more information see myFunction" help, publish MATLAB Answers — New Questions
Check if two colors are similar
How can I compare two RGB colors to see if they look the same (to human eyes). And I’d like to be able to control the amount of similarity.
Say for example 15. What is the best way to compare the colors to have more colors in the range.
I have two idea in my mind but I don’t know if they are the best options:
% (colorN(1) to colorN(3) are the RGB value of the color)
% Method 1:
if abs(double(color1(1)) – double(color2(1))) < 15 …
&& abs(double(color1(2)) – double(color2(2))) < 15 …
&& abs(double(color1(3)) – double(color2(3))) < 15
% The colors are similar
end
% Method 2
if sum(abs(double(color1(1)) – double(color2(1))), …
abs(double(color1(2)) – double(color2(2))), …
abs(double(color1(3)) – double(color2(3)))) < 15
% The colors are similar
end
% Maybe method 3 which is a combination of the two.
I’ve tried these a little bit by myself but I wasn’t happy with the results. Anyone knows a better way to compare? Does Matlab have a function for that?How can I compare two RGB colors to see if they look the same (to human eyes). And I’d like to be able to control the amount of similarity.
Say for example 15. What is the best way to compare the colors to have more colors in the range.
I have two idea in my mind but I don’t know if they are the best options:
% (colorN(1) to colorN(3) are the RGB value of the color)
% Method 1:
if abs(double(color1(1)) – double(color2(1))) < 15 …
&& abs(double(color1(2)) – double(color2(2))) < 15 …
&& abs(double(color1(3)) – double(color2(3))) < 15
% The colors are similar
end
% Method 2
if sum(abs(double(color1(1)) – double(color2(1))), …
abs(double(color1(2)) – double(color2(2))), …
abs(double(color1(3)) – double(color2(3)))) < 15
% The colors are similar
end
% Maybe method 3 which is a combination of the two.
I’ve tried these a little bit by myself but I wasn’t happy with the results. Anyone knows a better way to compare? Does Matlab have a function for that? How can I compare two RGB colors to see if they look the same (to human eyes). And I’d like to be able to control the amount of similarity.
Say for example 15. What is the best way to compare the colors to have more colors in the range.
I have two idea in my mind but I don’t know if they are the best options:
% (colorN(1) to colorN(3) are the RGB value of the color)
% Method 1:
if abs(double(color1(1)) – double(color2(1))) < 15 …
&& abs(double(color1(2)) – double(color2(2))) < 15 …
&& abs(double(color1(3)) – double(color2(3))) < 15
% The colors are similar
end
% Method 2
if sum(abs(double(color1(1)) – double(color2(1))), …
abs(double(color1(2)) – double(color2(2))), …
abs(double(color1(3)) – double(color2(3)))) < 15
% The colors are similar
end
% Maybe method 3 which is a combination of the two.
I’ve tried these a little bit by myself but I wasn’t happy with the results. Anyone knows a better way to compare? Does Matlab have a function for that? color, image processing, image, compare, rgb MATLAB Answers — New Questions
What’s the best color model to extract texture features from?
Just like L*a*b is the best color model to extract color features from images because it matches human color perspective,what is the best color model to extract GLCM texture features from and whyJust like L*a*b is the best color model to extract color features from images because it matches human color perspective,what is the best color model to extract GLCM texture features from and why Just like L*a*b is the best color model to extract color features from images because it matches human color perspective,what is the best color model to extract GLCM texture features from and why color, texture, colormodel, image MATLAB Answers — New Questions
Teams Messaging Gets Autocorrect
Autocorrect Common Errors in Teams Chat and Channel Messages
After complaining about the recently announced price rises for Microsoft 365, it’s time to reflect that part of the reason why people buy Microsoft 365 is the expectation of new functionality. Although not guaranteed by Microsoft, it’s a truism that software doesn’t remain static for long. Bugs are fixed, but more importantly, new features are added to make the software more competitive.
At this point, Microsoft 365 really only competes with itself. By this I mean that the reason why Microsoft adds new functionality to its products is often to convince customers to upgrade to a higher-priced and more feature-rich license. For example, from Office 365 E3 ($26/month after the increase) to Microsoft 365 E5 ($60/month after the increase). Good logic can underpin the decision to upgrade. For example, a company that needs stronger compliance solutions to satisfy industry regulations might find that the additional Purview functionality licensed by Microsoft 365 E5 is exactly what they need.
And sometimes Microsoft upgrades the base software to satisfy customer requirements or simply fill in a gap that is so obvious that you wonder why it’s taken Microsoft so long to deliver a solution.
Teams Messaging Autocorrect Rolling out to Targeted Tenants
Take the case of message center notification MC1192251 (5 December 2025, Microsoft 365 roadmap item 534487) which proclaims that “Autocorrect is coming to Microsoft Teams compose. Commonly misspelled words will now be automatically corrected while composing messages.” The change applies to Teams desktop on Windows and MacOS and is rolling out to targeted release tenants now with the intention of reaching general availability in mid-to-late January 2026. GCC will get the change about a month later.
Sorry, did I just report that Teams will autocorrect spellings when composing messages? This is the communication vehicle that was going to take over from email that is only just going to autocorrect spellings in 2026, nine years after the product launch when the other Office applications have been able to autocorrect text for years. Yes, that’s exactly the situation. Once the update is distributed to clients, Teams will correct “commonly” misspelled words (meaning no slang, local argot, or technical terms that are not in common usage) as people compose chat or channel messages.
Spell Checking and Autocorrect
Teams has been able to spell check for years, including the ability to autodetect the language used in a message. Autocorrect means that the editor detects common mistyping errors, like “mantain” and automatically replaces the error with the correct value (maintain in this case). It’s still perfectly possible to include a bunch of misspellings in a message because Autocorrect doesn’t pick up every possible mistake. Outlook settings (Figure 1) allow you to configure Autocorrect with new values to check for (like changing Msoft to Microsoft, or correctly capitalizing SharePoint), but that facility isn’t in Teams.
Update: According to Microsoft support, Teams uses a “system-level Autocorrect dictionary rather than an app-specific one, which helps maintain consistency across devices and languages” and that’s why individual users can’t add custom Autocorrect entries. I think that’s a flawed decision because everyone probably has some words they mistype frequently that won’t appear in a system-level dictionary.
In the case of Teams, autocorrect is enabled for all users. If someone wants to turn autocorrect off, they go to the General section of the Settings app and disable the option to Correct words while typing (Figure 2). Notice that there’s no way to add customize misspellings.
Little Things Mean a Lot
It’s good that Autocorrect has turned up in Teams messaging, even if it would be much better if individual users could configure Autocorrect to check for the misspellings that everyone is prone to make when typing.
In passing, let me note that the option to create an audio-only recording of a Teams meeting is now rolling out in production. It’s another example of an update that might not mean much to some but is extraordinarily useful in certain circumstances. You just need to be in the right situation!
Learn how to exploit the data available to Microsoft 365 tenant administrators through the Office 365 for IT Pros eBook. We love figuring out how things work.
Hide (or remove) the geoaxes in geobasemap
How can I hide or remove the geoaxes in my plot? When I use the command:
geobasemap(‘grayland’)
everything works as expected. Indeed, both the map and the axes are visible. However, when I try to hide the axes using the following command, nothing changes:
geobasemap(‘grayland’)
gx = geoaxes;
set(gx,’Visible’,’off’)How can I hide or remove the geoaxes in my plot? When I use the command:
geobasemap(‘grayland’)
everything works as expected. Indeed, both the map and the axes are visible. However, when I try to hide the axes using the following command, nothing changes:
geobasemap(‘grayland’)
gx = geoaxes;
set(gx,’Visible’,’off’) How can I hide or remove the geoaxes in my plot? When I use the command:
geobasemap(‘grayland’)
everything works as expected. Indeed, both the map and the axes are visible. However, when I try to hide the axes using the following command, nothing changes:
geobasemap(‘grayland’)
gx = geoaxes;
set(gx,’Visible’,’off’) geobasemap, geoaxes MATLAB Answers — New Questions
PID control simulation code improvement
I’m trying to simulate the roll angle dynamics using the PID controller. I just would like to know if this code iw well-written or can be improved?
Any help will be appreciated.
%% Parameters
J = 0.03; % moment of inertia (kg*m^2)
Kp = 2;
Ki = 1;
Kd = 0.5;
dt = 0.01; % time step
T = 3; % total time
t = 0:dt:T;
N = length(t);
phi_ref_deg = 5;
phi_ref = deg2rad(phi_ref_deg);
phi = zeros(1, N);
omega = zeros(1, N);
u = zeros(1, N);
e = zeros(1, N);
I = 0;
prev_e = phi_ref – phi(1);
for k = 1:N-1
e(k) = phi_ref – phi(k);
I = I + e(k)*dt;
D = (e(k) – prev_e)/dt;
% PID control
u(k) = Kp*e(k) + Ki*I + Kd*D;
% Dynamics: J*phi_ddot = u
phi_ddot = u(k)/J;
omega(k+1) = omega(k) + phi_ddot*dt;
phi(k+1) = phi(k) + omega(k+1)*dt;
prev_e = e(k);
end
e(end) = phi_ref – phi(end);
phi_deg = rad2deg(phi);
peak_val = max(phi_deg);
overshoot = max(0, peak_val – phi_ref_deg);
fprintf(‘Overshoot: %.3f degn’, overshoot);
figure;
subplot(2,1,1)
plot(t, phi_deg, ‘LineWidth’, 1.4); hold on;
yline(phi_ref_deg, ‘–r’, ‘Ref’);
xlabel(‘Time (s)’);
ylabel(‘Roll angle (deg)’);
title(‘Roll Angle Response’);
grid on;
subplot(2,1,2)
plot(t, u, ‘LineWidth’, 1.4);
xlabel(‘Time (s)’);
ylabel(‘u (N*m)’);
title(‘Control Torque’);
grid on;I’m trying to simulate the roll angle dynamics using the PID controller. I just would like to know if this code iw well-written or can be improved?
Any help will be appreciated.
%% Parameters
J = 0.03; % moment of inertia (kg*m^2)
Kp = 2;
Ki = 1;
Kd = 0.5;
dt = 0.01; % time step
T = 3; % total time
t = 0:dt:T;
N = length(t);
phi_ref_deg = 5;
phi_ref = deg2rad(phi_ref_deg);
phi = zeros(1, N);
omega = zeros(1, N);
u = zeros(1, N);
e = zeros(1, N);
I = 0;
prev_e = phi_ref – phi(1);
for k = 1:N-1
e(k) = phi_ref – phi(k);
I = I + e(k)*dt;
D = (e(k) – prev_e)/dt;
% PID control
u(k) = Kp*e(k) + Ki*I + Kd*D;
% Dynamics: J*phi_ddot = u
phi_ddot = u(k)/J;
omega(k+1) = omega(k) + phi_ddot*dt;
phi(k+1) = phi(k) + omega(k+1)*dt;
prev_e = e(k);
end
e(end) = phi_ref – phi(end);
phi_deg = rad2deg(phi);
peak_val = max(phi_deg);
overshoot = max(0, peak_val – phi_ref_deg);
fprintf(‘Overshoot: %.3f degn’, overshoot);
figure;
subplot(2,1,1)
plot(t, phi_deg, ‘LineWidth’, 1.4); hold on;
yline(phi_ref_deg, ‘–r’, ‘Ref’);
xlabel(‘Time (s)’);
ylabel(‘Roll angle (deg)’);
title(‘Roll Angle Response’);
grid on;
subplot(2,1,2)
plot(t, u, ‘LineWidth’, 1.4);
xlabel(‘Time (s)’);
ylabel(‘u (N*m)’);
title(‘Control Torque’);
grid on; I’m trying to simulate the roll angle dynamics using the PID controller. I just would like to know if this code iw well-written or can be improved?
Any help will be appreciated.
%% Parameters
J = 0.03; % moment of inertia (kg*m^2)
Kp = 2;
Ki = 1;
Kd = 0.5;
dt = 0.01; % time step
T = 3; % total time
t = 0:dt:T;
N = length(t);
phi_ref_deg = 5;
phi_ref = deg2rad(phi_ref_deg);
phi = zeros(1, N);
omega = zeros(1, N);
u = zeros(1, N);
e = zeros(1, N);
I = 0;
prev_e = phi_ref – phi(1);
for k = 1:N-1
e(k) = phi_ref – phi(k);
I = I + e(k)*dt;
D = (e(k) – prev_e)/dt;
% PID control
u(k) = Kp*e(k) + Ki*I + Kd*D;
% Dynamics: J*phi_ddot = u
phi_ddot = u(k)/J;
omega(k+1) = omega(k) + phi_ddot*dt;
phi(k+1) = phi(k) + omega(k+1)*dt;
prev_e = e(k);
end
e(end) = phi_ref – phi(end);
phi_deg = rad2deg(phi);
peak_val = max(phi_deg);
overshoot = max(0, peak_val – phi_ref_deg);
fprintf(‘Overshoot: %.3f degn’, overshoot);
figure;
subplot(2,1,1)
plot(t, phi_deg, ‘LineWidth’, 1.4); hold on;
yline(phi_ref_deg, ‘–r’, ‘Ref’);
xlabel(‘Time (s)’);
ylabel(‘Roll angle (deg)’);
title(‘Roll Angle Response’);
grid on;
subplot(2,1,2)
plot(t, u, ‘LineWidth’, 1.4);
xlabel(‘Time (s)’);
ylabel(‘u (N*m)’);
title(‘Control Torque’);
grid on; code improvement MATLAB Answers — New Questions
Open Excel file read name of data file, process data file write to six different excel output data files, Repeat in a for loop
for i=1:30
open excel INPUT file
read data fle names from Excel INPUT file
% internal loop
for j=1:6
process data
open one excel OUTPUT file
write results to excel file as a cell array along one row of excel output file
save excel file
close excel file
end
end
% have tried multiple methods, all of which failed, maybe not appropriate to MATLAB R2020b?for i=1:30
open excel INPUT file
read data fle names from Excel INPUT file
% internal loop
for j=1:6
process data
open one excel OUTPUT file
write results to excel file as a cell array along one row of excel output file
save excel file
close excel file
end
end
% have tried multiple methods, all of which failed, maybe not appropriate to MATLAB R2020b? for i=1:30
open excel INPUT file
read data fle names from Excel INPUT file
% internal loop
for j=1:6
process data
open one excel OUTPUT file
write results to excel file as a cell array along one row of excel output file
save excel file
close excel file
end
end
% have tried multiple methods, all of which failed, maybe not appropriate to MATLAB R2020b? open read write to excel file in loop MATLAB Answers — New Questions
Why does MATLAB get stuck in the “Initializing” or “Busy” state or take a long time to start?
Not so much a question as as add to another post (that I dont have enough reputation points to post to). If someone can add this to that post, I feel like it would be helpful.
This post was helpfull to me, untill last time it didn’t fix the issue. https://www.mathworks.com/matlabcentral/answers/92566-why-does-matlab-get-stuck-in-the-initializing-or-busy-state-or-take-a-long-time-to-start
After reaching out to Matlab Support I was able to get Matalb working again. This is the fix they provided me:
Thank you for contacting MathWorks Support.
This could be an uncommon issue with MATLAB connecting with the default printer driver, which can cause MATLAB to be stuck in Initializing for an extended period of time.
Slow initialization between "init desktop" and "psParser" phases has been known to be caused by a slow connection between MATLAB and the default printer. If this is the root cause, then the following workarounds are available:
A) Set your default printer to "Print to PDF", then see if the issue is resolved.
The link below is a guide from Microsoft to change the default printer.
https://support.microsoft.com/en-us/windows/set-a-default-printer-in-windows-e10cf8b8-e596-b102-bf84-c41022b5036f#ID0EBF=Windows_10
B) If you have printers on your network that you are connected to, check to see if temporarily removing or disabling these devices helps.
Hope this helps someone.Not so much a question as as add to another post (that I dont have enough reputation points to post to). If someone can add this to that post, I feel like it would be helpful.
This post was helpfull to me, untill last time it didn’t fix the issue. https://www.mathworks.com/matlabcentral/answers/92566-why-does-matlab-get-stuck-in-the-initializing-or-busy-state-or-take-a-long-time-to-start
After reaching out to Matlab Support I was able to get Matalb working again. This is the fix they provided me:
Thank you for contacting MathWorks Support.
This could be an uncommon issue with MATLAB connecting with the default printer driver, which can cause MATLAB to be stuck in Initializing for an extended period of time.
Slow initialization between "init desktop" and "psParser" phases has been known to be caused by a slow connection between MATLAB and the default printer. If this is the root cause, then the following workarounds are available:
A) Set your default printer to "Print to PDF", then see if the issue is resolved.
The link below is a guide from Microsoft to change the default printer.
https://support.microsoft.com/en-us/windows/set-a-default-printer-in-windows-e10cf8b8-e596-b102-bf84-c41022b5036f#ID0EBF=Windows_10
B) If you have printers on your network that you are connected to, check to see if temporarily removing or disabling these devices helps.
Hope this helps someone. Not so much a question as as add to another post (that I dont have enough reputation points to post to). If someone can add this to that post, I feel like it would be helpful.
This post was helpfull to me, untill last time it didn’t fix the issue. https://www.mathworks.com/matlabcentral/answers/92566-why-does-matlab-get-stuck-in-the-initializing-or-busy-state-or-take-a-long-time-to-start
After reaching out to Matlab Support I was able to get Matalb working again. This is the fix they provided me:
Thank you for contacting MathWorks Support.
This could be an uncommon issue with MATLAB connecting with the default printer driver, which can cause MATLAB to be stuck in Initializing for an extended period of time.
Slow initialization between "init desktop" and "psParser" phases has been known to be caused by a slow connection between MATLAB and the default printer. If this is the root cause, then the following workarounds are available:
A) Set your default printer to "Print to PDF", then see if the issue is resolved.
The link below is a guide from Microsoft to change the default printer.
https://support.microsoft.com/en-us/windows/set-a-default-printer-in-windows-e10cf8b8-e596-b102-bf84-c41022b5036f#ID0EBF=Windows_10
B) If you have printers on your network that you are connected to, check to see if temporarily removing or disabling these devices helps.
Hope this helps someone. initializing, busy, startup, printer communication problem MATLAB Answers — New Questions
numerical Instabilities for bessel functions
I write the following code but i used the following parametrs and i did not get any numerical instabilities:
RI = 34.5; % Inner radius (m)
ratio = 47.5/34.5; % Ratio R_E / R_I
RE = ratio*RI; % outer radius (m)
h = 200/34.5 * RI; % Water depth (m)
d = 38/34.5 * RI; % Interface depth (m)
b = h-d;
But when i switch the parameters with this then i get numerical instabilities
RI = 0.4; % Inner radius (m)
RE = 0.8; % Outer radius (m)
d = 0.3; % Draft (m)
h = 1.0; % Water depth (m)
b = h-d; %from cylinder bottom to seabed (m)
the following is code
clear all;
close all;
clc;
tic;
% diary(‘iter_log.txt’);
% diary on
% fprintf(‘Run started: %sn’, datestr(now));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters from your setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 10; % Number of N modes
M = 10; % Number of M modes
Q = 10; % Number of Q modes
RI = 0.4; % Inner radius (m)
RE = 0.8; % Outer radius (m)
d = 0.3; % Draft (m)
h = 1.0; % Water depth (m) % Interface depth (m)
b = h-d; % Distance from cylinder bottom to seabed (m)
g = 9.81; % Gravity (m/s^2)
tau =0.2; % porosity ratio
gamma = 1.0; % discharge coefficient
b1 = (1 – tau) / (2 * gamma * tau^2); % nonlinear coeff
tol = 1e-4; % convergence tolerance
max_iter = 100; % max iterations
l = 0; % Azimuthal order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === Compute Bessel roots and scaled chi values ===
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
roots = bessel0j(l, Q); % Zeros of J_l(x)
chi = roots ./ RI; % chi_q^l = roots scaled by radius
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ===== Define Range for k_0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X_kRE = linspace(0.0025, 0.16,100); %%%%%% this is now k_0
X_kRE = linspace(0.05, 4.5, 40);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%===============Initialize omega array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta0 = zeros(size(X_kRE)); % Preallocate
omega = zeros(size(X_kRE)); %%%% omega preallocate
% iters_used = zeros(length(X_kRE),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for idx = 1:length(X_kRE)
wk = X_kRE(idx); % dimensional wavenumber
omega(idx) = sqrt(g * wk * tanh(wk * h)); % dispersion relation
% fprintf(‘n— idx %3d (T=%6.3f s) —n’, idx, 2*pi/omega(idx));
% drawnow;
%%%%%%%%%%%%% Compute group velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cg(idx) = (g * tanh(wk * h) + g * wk * h * (sech(wk * h))^2) …
* omega(idx) / (2 * g * wk * tanh(wk * h));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% =======Compute a1 based on tau
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a1(idx) = 0.93 * (1 – tau) * Cg(idx);
% dissip(idx) = a1(idx)/omega(idx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%=============== Find derivative of Z0 at z = -d
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z0d = (cosh(wk * h) * wk * sinh(wk * b)) / (2 * wk * h + sinh(2 * wk * h));
fun_alpha = @(x) omega(idx)^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes
opts_lsq = optimoptions(‘lsqnonlin’,’Display’,’off’); % ← add this
for n = 1:N
if n == 1
k(n) = -1i * wk;
else
x0_left = (2*n – 3) * pi / (2*h);
x0_right = (2*n – 1) * pi / (2*h);
guess = (x0_left + x0_right)/2;
% Use lsqnonlin for better convergence
k(n) = lsqnonlin(fun_alpha, guess, x0_left, x0_right,opts_lsq);
%%%%%%%%%%%%% Derivative of Z_n at z = -d %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Znd(n) = -k(n) * (cos(k(n)*h) * sin(k(n)*b)) / (2*k(n)*h + sin(2*k(n)*h));
end
end
%% finding the root lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;
end
% %%%%%%%%%%%%%%%%%%%%% Define matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = zeros(M,N);
A(1,1) = (cosh(wk*h)*sinh(wk*b))/(wk*b*(2*wk*h+sinh(2*wk*h)))*(besselh(l, wk*RE)) /(dbesselh(l, wk*RE)); % Set A_00
for j = 2:N
A(1,j) =(cos(k(j)*h)*sin(k(j)*b) )/ (k(j)*b*(2*k(j)*h + sin(2*k(j)*h)))*( besselk(l, k(j)*RE) )…
/(dbesselk(l, k(j)*RE)); % Set A_0n
end
for i = 2:M
A(i,1) =(2 * cosh(wk*h) * (-1)^(i-1) * wk*sinh(wk*b)) / (b * (2*wk*h + sinh(2*wk*h))*(wk^2 + m(i)^2))…
* (besselh(l, wk*RE))/ (dbesselh(l, wk*RE)); % Set A_m0
end
for i = 2:M
for j = 2:N
A(i,j) = (2*cos(k(j)*h)*(-1)^(i-1) *k(j)* sin(k(j)*b))/ (b * (2*k(j)*h + sin(2*k(j)*h))*(k(j)^2-m(i)^2))…
*(besselk(l, k(j)*RE)) / (dbesselk(l, k(j)*RE));% Set A_mn
end
end
%%%%%%%%%%%%%%%%%%%%%%%% DefineMatrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N,M);
B(1,1) = (4*sinh(wk*b)) / (RE * wk*log(RE/RI) *cosh(wk*h)); %set B_00
%%%%%%%%%%%B_0m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(1,j) = Rml_prime_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%%%%%%%%B_n0%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
B(i,1)=(4*sin(k(i)*b))/(RE*k(i)*log(RE/RI)*cos(k(i)*h));% Set B_n0
end
%%%%%%%%%%%%%%%%%%%%%%%%%%B_nm%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
for j=2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(i,j) = Rml_prime_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * wk*log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for j = 2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(1,j) = Rml_prime_star_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
% %%%%%%% Define Cn0%%%%%%
for i = 2:N
C(i,1) = -4 * sin(k(i)*b) / (RE *k(i)* log(RE/RI) * cos(k(i)*h));
end
% % %%%%%% Define Cnm%%%%%%
for i = 2:N
for j =2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(i,j) = Rml_prime_star_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%% write Matrix D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = zeros(M,N);
%%% write D00 %%%%%%
D(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2*wk*h + sinh(2*wk*h))) * (besselj(l, wk*RI))/ (dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for j =2:N
D(1,j) = (cos(k(j)*h) * sin(k(j)*b)) / (k(j)*b * (2*k(j)*h + sin(2*k(j)*h))) * (besseli(l, k(j)*RI) )/ (dbesseli(l, k(j)*RI));
end
%%%%%% write Dm0 %%%%%%%%%
for i = 2:M
D(i,1) = (2 * cosh(wk*h) * (-1)^(i-1) * wk * sinh(wk*b)) /(b * (2*wk*h + sinh(2*wk*h)) * (wk^2 + m(i)^2))…
*(besselj(l, wk*RI) )/(dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for i = 2:M
for j = 2:N
D(i,j) = (2 * cos(k(j)*h) * (-1)^(i-1) * k(j) * sin(k(j)*b)) /(b * (2*k(j)*h + sin(2*k(j)*h)) * (k(j)^2 – m(i)^2))…
*(besseli(l, k(j)*RI)) / (dbesseli(l, k(j)*RI));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M); % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
E(1,1) = (4 * sinh(wk*b)) / (RI *wk* log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(1,j) = Rml_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
%
end
%%%%%% Define Eno%%%%%%%%%%%%%%
for i = 2:N
E(i,1) = (4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
for i = 2:N
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(i,j) = Rml_prime_RI * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * wk*log(RE/RI) * cosh(wk*h));
% %%%%%% Define F0m%%%%
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(1,j) = Rml_star_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%% Defin Fn0%%%%%%
for i = 2:N
F(i,1) = (-4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
%%%%%%% Define Fnm %%%%%%%
for i = 2:N
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(i,j) = Rml_star_prime_RI * (4 * k(i)* (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
% % %%%%%%%%%% Write Matrix W %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = zeros(N,Q);
% %%W_{0q}
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 – r); % coth(chi*b) via decaying exp
W(1,q) = -(4*chi(q)/cosh(wk*h))*(( wk*sinh(wk*b)*coth_chib – chi(q)*cosh(wk*b) ) / ( wk^2 – chi(q)^2 ) ); %prof.chen
%
% % W(1, q) = – (4 * chi(q) / cosh(wk * h)) * (chi(q) * sinh(chi(q) * b) * cosh(wk * b) – wk * sinh(wk * b) * cosh(chi(q) * b)) …
% % / ((chi(q)^2 – wk^2) * sinh(chi(q) * b));
end
% % %%% Write W_{nq}
for n = 2:N
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 – r); % coth(chi*b), stable
W(n,q) = -(4*chi(q)/cos(k(n)*h)) * …
( ( k(n)*sin(k(n)*b)*coth_chib + chi(q)*cos(k(n)*b) ) … %prof.chen
/ ( k(n)^2 + chi(q)^2 ) );
% % W(n, q) = – (4 * chi(q) / cos(k(n) * h)) * (chi(q) * sinh(chi(q) * b) * cos(k(n) * b) + k(n) * sin(k(n) * b) * cosh(chi(q) * b)) …
% % / ((chi(q)^2 + k(n)^2) * sinh(chi(q) * b));
end
end
%%%%%%%%%%%%%%%%%%%%%%%Find E1_{q}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f2 = omega(idx)^2 / g;
E1 = zeros(Q,1);
S_q = zeros(Q,1);
for q = 1:Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% this expression is same as prof. Chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
exp2chi_d = exp(-2 * chi(q) * d); % e^{-2*chi*d}
exp2chi_b = exp(-2 * chi(q) * b); % e^{-2*chi*(h-d)}
exp2chi_h = exp(-2 * chi(q) * h); % e^{-2*chi*h}
exp1chi_d = exp(-chi(q) * d); % e^{-chi*d}
% intermediate ratio (same order as Fortran)
BB = (f2 – chi(q)) * (exp2chi_b / exp2chi_d);
BB = BB + (f2 + chi(q)) * (1 + 2 * exp2chi_b);
BB = BB / (f2 – chi(q) + (f2 + chi(q)) * exp2chi_d);
BB = BB / (1 + exp2chi_b);
% final coefficient
E1(q) = 2 * (BB * exp2chi_d – 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This expression is that one derived directly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% % N0, D0, H0 exactly as in the math above
% N0 = (chi(q) – f2) + (chi(q) + f2) * exp2chi_h;
% D0 = (f2 – chi(q)) + (f2 + chi(q)) * exp2chi_d;
% H0 = 1 + exp2chi_b;
%
% % final coefficient (no add/subtract step)
% E1(q) = 2 * N0 / (D0 * H0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This one is the hyperbolic one form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% num = chi(q)*cosh(chi(q)*h) – f2*sinh(chi(q)*h);
% den = (f2*cosh(chi(q)*d) – chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*(h – d)); % or: * cosh(chi(q)*b)
% E1(q) = num / den;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% This is the upper region velocity potential at z=0, write
%%%%%%%%%%%%%%%% from prof.chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CCmc = 1 + (exp2chi_b / exp2chi_d) + 2*exp2chi_b;
CCmc = CCmc / (1 + exp2chi_b);
CCmc = BB * (1 + exp2chi_d) – CCmc;
S_q(q) = CCmc * exp1chi_d; % this equals the summand value for this q
% — direct hyperbolic (exact) —
% % % C_hyp = ( chi(q)*cosh(chi(q)*h) – f2*sinh(chi(q)*h) ) …
% % % / ( (f2*cosh(chi(q)*d) – chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*b) );
% % %
% % % S_hyp(q) = C_hyp * cosh(chi(q)*d) + sinh(chi(q)*h) / cosh(chi(q)*b); % summand
end
% H(q): diagonal term (exponential form)
H1 = zeros(Q,1);
for q = 1:Q
% dissip =
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
mid = 4*r/(1 – r^2) – E1(q); % 2/sinh(2*chi*b) – C_m
H1(q) = mid + 1i*a1(idx)*chi(q)/omega(idx); % mid + i*(a1/w)*chi
end
H = diag(H1);
%%%%%%%%%%matric G%%%%%%%%%%%%%%%%%%
G = zeros(Q, N); % Preallocate
for q = 1:Q
G(q,1) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Z0d …
* ( chi(q) / (chi(q)^2 – wk^2) )* ( besselj(l, wk*RI) / dbesselj(l, wk*RI) ); % this one is my calculated
end
% %G_qn
for q = 1:Q
for n = 2:N
G(q, n) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Znd(n) …
* ( chi(q) / (k(n)^2 + chi(q)^2) )* ( besseli(l, k(n)*RI) / dbesseli(l, k(n)*RI) ); % this one is my calculated
end
end
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Define the right hand side vector
U = zeros(2*M + 2*N + Q, 1); % Full vector, M+N+M+N+Q elements
% % Block 1: U (size M = 4)
for i = 1:M
if i == 1
U(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
U(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2)*cosh(wk*h))*besselj(l, wk*RE); %Z_m^l
end
end
% % Block 2: Y (size N = 3)
for j = 1:N
if j == 1
U(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
U(j + M, 1) = 0; % Y_n^l
end
end
% % Block 3: X (size M = 4)
for i = 1:M
U(i + M + N, 1) = 0; % X_0^l, X_m^l
end
% Block 4: W (size N = 3)
for j = 1:N
U(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
% % Identity matrices
I_M = eye(M);
I_N = eye(N);
% Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);
ZMQ = zeros(M,Q);
ZNQ = zeros(N,Q);
ZQM = zeros(Q,M);
ZQN = zeros(Q,N);
% Construct the full block matrix S
S = [ I_M, -A, ZMM, ZMN, ZMQ;
-B, I_N, -C, ZNN, ZNQ;
ZMM, ZMN, I_M, -D, ZMQ;
-E, ZNN, -F, I_N, -W;
ZQM, ZQN, ZQM, -G, H ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
converged = false;
psi = zeros(Q,1);
T_old = []; % <– track the full solution from prev. iter
% T_old = zeros(2*M + 2*N + Q, 1); % previous unknowns (seed)
for iter = 1:max_iter
% (A) overwrite ONLY Block 5 of U from current psi
for q = 1:Q
U(2*M + 2*N + q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi(q);
end
T = S U; % Solve the linear system
% T = pinv(S) * U; % Use pseudoinverse for stability
b_vec = T(1:M); % Coefficients b^l
a_vec = T(M+1 : M+N); % Coefficients a^l
c_vec = T(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = T(2*M+N+1:2*M+2*N); % Coefficients d^l
e_vec = T(2*M+2*N+1:end); % (Q×1)
% (D) update psi for NEXT iteration from CURRENT coefficients
for q = 1:Q
integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) …
.* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) …
.* besselj(l, chi(q)*r) .* r;
psi(q) = integral(integrand, 0, RI, ‘AbsTol’,1e-8, ‘RelTol’,1e-6);
end
% % % drawnow; % optional: flush output each iter
% === compact per-iteration print, only for T in [5,35] ===
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf(‘idx %3d | iter %2d | ||psi|| = %.3en’, idx, iter, norm(psi));
% drawnow;
% end
% % % % (E) convergence on ALL unknowns (T vs previous T)
% — per-iteration diagnostics (optional) —
if ~isempty(T_old)
diff_T = max(abs(T – T_old));
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf(‘k0 idx %3d | iter %2d: max|ΔT| = %.3e, ||psi|| = %.3en’, …
% idx, iter, diff_T, norm(psi));
% drawnow; % ← add this to flush each iteration line too
% end
if diff_T < tol
converged = true;
break
end
end
T_old = T;
end % <<< end of iter loop
% % —- ONE summary line per frequency (place it here) —-
% iters_used(idx) = iter; % how many iterations this frequency needed
% fprintf(‘idx %3d (T=%6.3f s): iters=%2dn’, idx, 2*pi/omega(idx), iter);
% drawnow;
% ————– end nonlinear iteration ————–
% % %%%%%%%%% Wave motion%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
term1 = (d_vec(1)*cosh(wk*h)^2) / ((2*wk*h + sinh(2*wk*h)) * dbesselj(l, wk*RI));
sum1 = 0;
for n =2:N
sum1 = sum1 + d_vec(n)*cos(k(n)*h)^2/ ((2*k(n)*h + sin(2*k(n)*h))* dbesseli(l, k(n)*RI));
end
phiUell = 0;
for q =1:Q
phiUell = phiUell +e_vec(q)* S_q(q);
end
eta0(idx) = abs(term1+sum1+phiUell);
end
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % Plotting Data
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% figure(2); clf;
% >>> CHANGES START (3): one figure = two taus + Dr. Chen only
% plot(2*pi./omega, eta0, ‘k’, ‘LineWidth’, 1.5); %%% this is T which is T= 2*pi/omega plotting against T
plot(omega.^2 * RE / g, eta0, ‘k’, ‘LineWidth’, 1.5); %%% this is T which is T= 2*pi/omega plotting against T
hold on; % Add extra plots on same figure
% % % %
% % % % % Now plot the 3 CSV experimental points
% scatter(data_chen(:,1), data_chen(:,2), 30, ‘r’, ‘o’, ‘filled’);
% % scatter(data_liu(:,1), data_liu(:,2), 30, ‘g’, ‘s’, ‘filled’);
% scatter(data_exp(:,1), data_exp(:,2), 30, ‘b’, ‘^’, ‘filled’);
xlabel(‘$T$’, ‘Interpreter’, ‘latex’);
ylabel(‘$|eta / (iA)|$’, ‘Interpreter’, ‘latex’);
title(‘Wave motion amplitude for $R_E = 138$’, ‘Interpreter’, ‘latex’);
legend({‘$tau=0.2$’,’Model test’}, …
‘Interpreter’,’latex’,’Location’,’northwest’);
grid on;
% xlim([5 35]); % Match the reference plot
% ylim([0 4.5]); % Optional: based on expected peak height
xlim([0 4]); % Match the reference plot
ylim([0 7]); % Optional: based on expected peak height
% === plotting section ===
% diary off
elapsedTime = toc;
disp([‘Time consuming = ‘, num2str(elapsedTime), ‘ s’]);
function out = dbesselk(l, z)
%DBESSELK Derivative of the modified Bessel function of the second kind
% out = dbesselk(l, z)
% Returns d/dz [K_l(z)] using the recurrence formula:
% K_l'(z) = -0.5 * (K_{l-1}(z) + K_{l+1}(z))
out = -0.5 * (besselk(l-1, z) + besselk(l+1, z));
end
function out = dbesselj(l, z)
%DBESSELJ Derivative of the Bessel function of the first kind
% out = dbesselj(l, z)
% Returns d/dz [J_l(z)] using the recurrence formula:
% J_l'(z) = 0.5 * (J_{l-1}(z) – J_{l+1}(z))
out = 0.5 * (besselj(l-1, z) – besselj(l+1, z));
end
function out = dbesseli(l, z)
%DBESSELI Derivative of the modified Bessel function of the first kind
% out = dbesseli(l, z)
% Returns d/dz [I_l(z)] using the recurrence formula:
% I_l'(z) = 0.5 * (I_{l-1}(z) + I_{l+1}(z))
out = 0.5 * (besseli(l-1, z) + besseli(l+1, z));
end
function out = dbesselh(l, z)
%DBESSELH Derivative of the Hankel function of the first kind
% out = dbesselh(l, z)
% Returns d/dz [H_l^{(1)}(z)] using the recurrence formula:
% H_l^{(1)’}(z) = 0.5 * (H_{l-1}^{(1)}(z) – H_{l+1}^{(1)}(z))
out = 0.5 * (besselh(l-1, 1, z) – besselh(l+1, 1, z));
end
function x = bessel0j(l,q,opt)
% a row vector of the first q roots of bessel function Jl(x), integer order.
% if opt = ‘d’, row vector of the first q roots of dJl(x)/dx, integer order.
% if opt is not provided, the default is zeros of Jl(x).
% all roots are positive, except when l=0,
% x=0 is included as a root of dJ0(x)/dx (standard convention).
%
% starting point for for zeros of Jl was borrowed from Cleve Moler,
% but the starting points for both Jl and Jl’ can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13.
%
% David Goodmanson
%
% x = bessel0j(l,q,opt)
k = 1:q;
if nargin==3 && opt==’d’
beta = (k + l/2 – 3/4)*pi;
mu = 4*l^2;
x = beta – (mu+3)./(8*beta) – 4*(7*mu^2+82*mu-9)./(3*(8*beta).^3);
for j=1:8
xnew = x – besseljd(l,x)./ …
(besselj(l,x).*((l^2./x.^2)-1) -besseljd(l,x)./x);
x = xnew;
end
if l==0
x(1) = 0; % correct a small numerical difference from 0
end
else
beta = (k + l/2 – 1/4)*pi;
mu = 4*l^2;
x = beta – (mu-1)./(8*beta) – 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);
for j=1:8
xnew = x – besselj(l,x)./besseljd(l,x);
x = xnew;
end
end
end
% — Local helper function for derivative of Bessel function —
function dJ = besseljd(l, x)
dJ = 0.5 * (besselj(l – 1, x) – besselj(l + 1, x));
end
function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi)
% n = 0 mode
term1 = (d_vec(1) * Z0d * besselj(l, wk*r)) / (dbesselj(l, wk*RI));
% n >= 2 evanescent sum
sum1 = 0;
for nidx = 2:N
sum1 = sum1 + d_vec(nidx) * Znd(nidx) * besseli(l, k(nidx)*r) …
/( dbesseli(l, k(nidx)*RI));
end
% q = 1..Q radial sum
sum2 = 0;
for qidx = 1:Q
sum2 = sum2 + e_vec(qidx) *chi(qidx)* besselj(l, chi(qidx)*r) / dbesselj(l, chi(qidx)*RI);
end
v_D_val = term1 + sum1 + sum2;
endI write the following code but i used the following parametrs and i did not get any numerical instabilities:
RI = 34.5; % Inner radius (m)
ratio = 47.5/34.5; % Ratio R_E / R_I
RE = ratio*RI; % outer radius (m)
h = 200/34.5 * RI; % Water depth (m)
d = 38/34.5 * RI; % Interface depth (m)
b = h-d;
But when i switch the parameters with this then i get numerical instabilities
RI = 0.4; % Inner radius (m)
RE = 0.8; % Outer radius (m)
d = 0.3; % Draft (m)
h = 1.0; % Water depth (m)
b = h-d; %from cylinder bottom to seabed (m)
the following is code
clear all;
close all;
clc;
tic;
% diary(‘iter_log.txt’);
% diary on
% fprintf(‘Run started: %sn’, datestr(now));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters from your setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 10; % Number of N modes
M = 10; % Number of M modes
Q = 10; % Number of Q modes
RI = 0.4; % Inner radius (m)
RE = 0.8; % Outer radius (m)
d = 0.3; % Draft (m)
h = 1.0; % Water depth (m) % Interface depth (m)
b = h-d; % Distance from cylinder bottom to seabed (m)
g = 9.81; % Gravity (m/s^2)
tau =0.2; % porosity ratio
gamma = 1.0; % discharge coefficient
b1 = (1 – tau) / (2 * gamma * tau^2); % nonlinear coeff
tol = 1e-4; % convergence tolerance
max_iter = 100; % max iterations
l = 0; % Azimuthal order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === Compute Bessel roots and scaled chi values ===
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
roots = bessel0j(l, Q); % Zeros of J_l(x)
chi = roots ./ RI; % chi_q^l = roots scaled by radius
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ===== Define Range for k_0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X_kRE = linspace(0.0025, 0.16,100); %%%%%% this is now k_0
X_kRE = linspace(0.05, 4.5, 40);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%===============Initialize omega array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta0 = zeros(size(X_kRE)); % Preallocate
omega = zeros(size(X_kRE)); %%%% omega preallocate
% iters_used = zeros(length(X_kRE),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for idx = 1:length(X_kRE)
wk = X_kRE(idx); % dimensional wavenumber
omega(idx) = sqrt(g * wk * tanh(wk * h)); % dispersion relation
% fprintf(‘n— idx %3d (T=%6.3f s) —n’, idx, 2*pi/omega(idx));
% drawnow;
%%%%%%%%%%%%% Compute group velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cg(idx) = (g * tanh(wk * h) + g * wk * h * (sech(wk * h))^2) …
* omega(idx) / (2 * g * wk * tanh(wk * h));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% =======Compute a1 based on tau
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a1(idx) = 0.93 * (1 – tau) * Cg(idx);
% dissip(idx) = a1(idx)/omega(idx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%=============== Find derivative of Z0 at z = -d
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z0d = (cosh(wk * h) * wk * sinh(wk * b)) / (2 * wk * h + sinh(2 * wk * h));
fun_alpha = @(x) omega(idx)^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes
opts_lsq = optimoptions(‘lsqnonlin’,’Display’,’off’); % ← add this
for n = 1:N
if n == 1
k(n) = -1i * wk;
else
x0_left = (2*n – 3) * pi / (2*h);
x0_right = (2*n – 1) * pi / (2*h);
guess = (x0_left + x0_right)/2;
% Use lsqnonlin for better convergence
k(n) = lsqnonlin(fun_alpha, guess, x0_left, x0_right,opts_lsq);
%%%%%%%%%%%%% Derivative of Z_n at z = -d %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Znd(n) = -k(n) * (cos(k(n)*h) * sin(k(n)*b)) / (2*k(n)*h + sin(2*k(n)*h));
end
end
%% finding the root lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;
end
% %%%%%%%%%%%%%%%%%%%%% Define matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = zeros(M,N);
A(1,1) = (cosh(wk*h)*sinh(wk*b))/(wk*b*(2*wk*h+sinh(2*wk*h)))*(besselh(l, wk*RE)) /(dbesselh(l, wk*RE)); % Set A_00
for j = 2:N
A(1,j) =(cos(k(j)*h)*sin(k(j)*b) )/ (k(j)*b*(2*k(j)*h + sin(2*k(j)*h)))*( besselk(l, k(j)*RE) )…
/(dbesselk(l, k(j)*RE)); % Set A_0n
end
for i = 2:M
A(i,1) =(2 * cosh(wk*h) * (-1)^(i-1) * wk*sinh(wk*b)) / (b * (2*wk*h + sinh(2*wk*h))*(wk^2 + m(i)^2))…
* (besselh(l, wk*RE))/ (dbesselh(l, wk*RE)); % Set A_m0
end
for i = 2:M
for j = 2:N
A(i,j) = (2*cos(k(j)*h)*(-1)^(i-1) *k(j)* sin(k(j)*b))/ (b * (2*k(j)*h + sin(2*k(j)*h))*(k(j)^2-m(i)^2))…
*(besselk(l, k(j)*RE)) / (dbesselk(l, k(j)*RE));% Set A_mn
end
end
%%%%%%%%%%%%%%%%%%%%%%%% DefineMatrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N,M);
B(1,1) = (4*sinh(wk*b)) / (RE * wk*log(RE/RI) *cosh(wk*h)); %set B_00
%%%%%%%%%%%B_0m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(1,j) = Rml_prime_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%%%%%%%%B_n0%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
B(i,1)=(4*sin(k(i)*b))/(RE*k(i)*log(RE/RI)*cos(k(i)*h));% Set B_n0
end
%%%%%%%%%%%%%%%%%%%%%%%%%%B_nm%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
for j=2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(i,j) = Rml_prime_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * wk*log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for j = 2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(1,j) = Rml_prime_star_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
% %%%%%%% Define Cn0%%%%%%
for i = 2:N
C(i,1) = -4 * sin(k(i)*b) / (RE *k(i)* log(RE/RI) * cos(k(i)*h));
end
% % %%%%%% Define Cnm%%%%%%
for i = 2:N
for j =2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(i,j) = Rml_prime_star_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%% write Matrix D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = zeros(M,N);
%%% write D00 %%%%%%
D(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2*wk*h + sinh(2*wk*h))) * (besselj(l, wk*RI))/ (dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for j =2:N
D(1,j) = (cos(k(j)*h) * sin(k(j)*b)) / (k(j)*b * (2*k(j)*h + sin(2*k(j)*h))) * (besseli(l, k(j)*RI) )/ (dbesseli(l, k(j)*RI));
end
%%%%%% write Dm0 %%%%%%%%%
for i = 2:M
D(i,1) = (2 * cosh(wk*h) * (-1)^(i-1) * wk * sinh(wk*b)) /(b * (2*wk*h + sinh(2*wk*h)) * (wk^2 + m(i)^2))…
*(besselj(l, wk*RI) )/(dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for i = 2:M
for j = 2:N
D(i,j) = (2 * cos(k(j)*h) * (-1)^(i-1) * k(j) * sin(k(j)*b)) /(b * (2*k(j)*h + sin(2*k(j)*h)) * (k(j)^2 – m(i)^2))…
*(besseli(l, k(j)*RI)) / (dbesseli(l, k(j)*RI));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M); % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
E(1,1) = (4 * sinh(wk*b)) / (RI *wk* log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(1,j) = Rml_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
%
end
%%%%%% Define Eno%%%%%%%%%%%%%%
for i = 2:N
E(i,1) = (4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
for i = 2:N
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(i,j) = Rml_prime_RI * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * wk*log(RE/RI) * cosh(wk*h));
% %%%%%% Define F0m%%%%
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(1,j) = Rml_star_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%% Defin Fn0%%%%%%
for i = 2:N
F(i,1) = (-4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
%%%%%%% Define Fnm %%%%%%%
for i = 2:N
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(i,j) = Rml_star_prime_RI * (4 * k(i)* (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
% % %%%%%%%%%% Write Matrix W %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = zeros(N,Q);
% %%W_{0q}
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 – r); % coth(chi*b) via decaying exp
W(1,q) = -(4*chi(q)/cosh(wk*h))*(( wk*sinh(wk*b)*coth_chib – chi(q)*cosh(wk*b) ) / ( wk^2 – chi(q)^2 ) ); %prof.chen
%
% % W(1, q) = – (4 * chi(q) / cosh(wk * h)) * (chi(q) * sinh(chi(q) * b) * cosh(wk * b) – wk * sinh(wk * b) * cosh(chi(q) * b)) …
% % / ((chi(q)^2 – wk^2) * sinh(chi(q) * b));
end
% % %%% Write W_{nq}
for n = 2:N
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 – r); % coth(chi*b), stable
W(n,q) = -(4*chi(q)/cos(k(n)*h)) * …
( ( k(n)*sin(k(n)*b)*coth_chib + chi(q)*cos(k(n)*b) ) … %prof.chen
/ ( k(n)^2 + chi(q)^2 ) );
% % W(n, q) = – (4 * chi(q) / cos(k(n) * h)) * (chi(q) * sinh(chi(q) * b) * cos(k(n) * b) + k(n) * sin(k(n) * b) * cosh(chi(q) * b)) …
% % / ((chi(q)^2 + k(n)^2) * sinh(chi(q) * b));
end
end
%%%%%%%%%%%%%%%%%%%%%%%Find E1_{q}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f2 = omega(idx)^2 / g;
E1 = zeros(Q,1);
S_q = zeros(Q,1);
for q = 1:Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% this expression is same as prof. Chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
exp2chi_d = exp(-2 * chi(q) * d); % e^{-2*chi*d}
exp2chi_b = exp(-2 * chi(q) * b); % e^{-2*chi*(h-d)}
exp2chi_h = exp(-2 * chi(q) * h); % e^{-2*chi*h}
exp1chi_d = exp(-chi(q) * d); % e^{-chi*d}
% intermediate ratio (same order as Fortran)
BB = (f2 – chi(q)) * (exp2chi_b / exp2chi_d);
BB = BB + (f2 + chi(q)) * (1 + 2 * exp2chi_b);
BB = BB / (f2 – chi(q) + (f2 + chi(q)) * exp2chi_d);
BB = BB / (1 + exp2chi_b);
% final coefficient
E1(q) = 2 * (BB * exp2chi_d – 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This expression is that one derived directly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% % N0, D0, H0 exactly as in the math above
% N0 = (chi(q) – f2) + (chi(q) + f2) * exp2chi_h;
% D0 = (f2 – chi(q)) + (f2 + chi(q)) * exp2chi_d;
% H0 = 1 + exp2chi_b;
%
% % final coefficient (no add/subtract step)
% E1(q) = 2 * N0 / (D0 * H0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This one is the hyperbolic one form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% num = chi(q)*cosh(chi(q)*h) – f2*sinh(chi(q)*h);
% den = (f2*cosh(chi(q)*d) – chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*(h – d)); % or: * cosh(chi(q)*b)
% E1(q) = num / den;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% This is the upper region velocity potential at z=0, write
%%%%%%%%%%%%%%%% from prof.chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CCmc = 1 + (exp2chi_b / exp2chi_d) + 2*exp2chi_b;
CCmc = CCmc / (1 + exp2chi_b);
CCmc = BB * (1 + exp2chi_d) – CCmc;
S_q(q) = CCmc * exp1chi_d; % this equals the summand value for this q
% — direct hyperbolic (exact) —
% % % C_hyp = ( chi(q)*cosh(chi(q)*h) – f2*sinh(chi(q)*h) ) …
% % % / ( (f2*cosh(chi(q)*d) – chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*b) );
% % %
% % % S_hyp(q) = C_hyp * cosh(chi(q)*d) + sinh(chi(q)*h) / cosh(chi(q)*b); % summand
end
% H(q): diagonal term (exponential form)
H1 = zeros(Q,1);
for q = 1:Q
% dissip =
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
mid = 4*r/(1 – r^2) – E1(q); % 2/sinh(2*chi*b) – C_m
H1(q) = mid + 1i*a1(idx)*chi(q)/omega(idx); % mid + i*(a1/w)*chi
end
H = diag(H1);
%%%%%%%%%%matric G%%%%%%%%%%%%%%%%%%
G = zeros(Q, N); % Preallocate
for q = 1:Q
G(q,1) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Z0d …
* ( chi(q) / (chi(q)^2 – wk^2) )* ( besselj(l, wk*RI) / dbesselj(l, wk*RI) ); % this one is my calculated
end
% %G_qn
for q = 1:Q
for n = 2:N
G(q, n) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Znd(n) …
* ( chi(q) / (k(n)^2 + chi(q)^2) )* ( besseli(l, k(n)*RI) / dbesseli(l, k(n)*RI) ); % this one is my calculated
end
end
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Define the right hand side vector
U = zeros(2*M + 2*N + Q, 1); % Full vector, M+N+M+N+Q elements
% % Block 1: U (size M = 4)
for i = 1:M
if i == 1
U(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
U(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2)*cosh(wk*h))*besselj(l, wk*RE); %Z_m^l
end
end
% % Block 2: Y (size N = 3)
for j = 1:N
if j == 1
U(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
U(j + M, 1) = 0; % Y_n^l
end
end
% % Block 3: X (size M = 4)
for i = 1:M
U(i + M + N, 1) = 0; % X_0^l, X_m^l
end
% Block 4: W (size N = 3)
for j = 1:N
U(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
% % Identity matrices
I_M = eye(M);
I_N = eye(N);
% Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);
ZMQ = zeros(M,Q);
ZNQ = zeros(N,Q);
ZQM = zeros(Q,M);
ZQN = zeros(Q,N);
% Construct the full block matrix S
S = [ I_M, -A, ZMM, ZMN, ZMQ;
-B, I_N, -C, ZNN, ZNQ;
ZMM, ZMN, I_M, -D, ZMQ;
-E, ZNN, -F, I_N, -W;
ZQM, ZQN, ZQM, -G, H ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
converged = false;
psi = zeros(Q,1);
T_old = []; % <– track the full solution from prev. iter
% T_old = zeros(2*M + 2*N + Q, 1); % previous unknowns (seed)
for iter = 1:max_iter
% (A) overwrite ONLY Block 5 of U from current psi
for q = 1:Q
U(2*M + 2*N + q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi(q);
end
T = S U; % Solve the linear system
% T = pinv(S) * U; % Use pseudoinverse for stability
b_vec = T(1:M); % Coefficients b^l
a_vec = T(M+1 : M+N); % Coefficients a^l
c_vec = T(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = T(2*M+N+1:2*M+2*N); % Coefficients d^l
e_vec = T(2*M+2*N+1:end); % (Q×1)
% (D) update psi for NEXT iteration from CURRENT coefficients
for q = 1:Q
integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) …
.* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) …
.* besselj(l, chi(q)*r) .* r;
psi(q) = integral(integrand, 0, RI, ‘AbsTol’,1e-8, ‘RelTol’,1e-6);
end
% % % drawnow; % optional: flush output each iter
% === compact per-iteration print, only for T in [5,35] ===
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf(‘idx %3d | iter %2d | ||psi|| = %.3en’, idx, iter, norm(psi));
% drawnow;
% end
% % % % (E) convergence on ALL unknowns (T vs previous T)
% — per-iteration diagnostics (optional) —
if ~isempty(T_old)
diff_T = max(abs(T – T_old));
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf(‘k0 idx %3d | iter %2d: max|ΔT| = %.3e, ||psi|| = %.3en’, …
% idx, iter, diff_T, norm(psi));
% drawnow; % ← add this to flush each iteration line too
% end
if diff_T < tol
converged = true;
break
end
end
T_old = T;
end % <<< end of iter loop
% % —- ONE summary line per frequency (place it here) —-
% iters_used(idx) = iter; % how many iterations this frequency needed
% fprintf(‘idx %3d (T=%6.3f s): iters=%2dn’, idx, 2*pi/omega(idx), iter);
% drawnow;
% ————– end nonlinear iteration ————–
% % %%%%%%%%% Wave motion%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
term1 = (d_vec(1)*cosh(wk*h)^2) / ((2*wk*h + sinh(2*wk*h)) * dbesselj(l, wk*RI));
sum1 = 0;
for n =2:N
sum1 = sum1 + d_vec(n)*cos(k(n)*h)^2/ ((2*k(n)*h + sin(2*k(n)*h))* dbesseli(l, k(n)*RI));
end
phiUell = 0;
for q =1:Q
phiUell = phiUell +e_vec(q)* S_q(q);
end
eta0(idx) = abs(term1+sum1+phiUell);
end
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % Plotting Data
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% figure(2); clf;
% >>> CHANGES START (3): one figure = two taus + Dr. Chen only
% plot(2*pi./omega, eta0, ‘k’, ‘LineWidth’, 1.5); %%% this is T which is T= 2*pi/omega plotting against T
plot(omega.^2 * RE / g, eta0, ‘k’, ‘LineWidth’, 1.5); %%% this is T which is T= 2*pi/omega plotting against T
hold on; % Add extra plots on same figure
% % % %
% % % % % Now plot the 3 CSV experimental points
% scatter(data_chen(:,1), data_chen(:,2), 30, ‘r’, ‘o’, ‘filled’);
% % scatter(data_liu(:,1), data_liu(:,2), 30, ‘g’, ‘s’, ‘filled’);
% scatter(data_exp(:,1), data_exp(:,2), 30, ‘b’, ‘^’, ‘filled’);
xlabel(‘$T$’, ‘Interpreter’, ‘latex’);
ylabel(‘$|eta / (iA)|$’, ‘Interpreter’, ‘latex’);
title(‘Wave motion amplitude for $R_E = 138$’, ‘Interpreter’, ‘latex’);
legend({‘$tau=0.2$’,’Model test’}, …
‘Interpreter’,’latex’,’Location’,’northwest’);
grid on;
% xlim([5 35]); % Match the reference plot
% ylim([0 4.5]); % Optional: based on expected peak height
xlim([0 4]); % Match the reference plot
ylim([0 7]); % Optional: based on expected peak height
% === plotting section ===
% diary off
elapsedTime = toc;
disp([‘Time consuming = ‘, num2str(elapsedTime), ‘ s’]);
function out = dbesselk(l, z)
%DBESSELK Derivative of the modified Bessel function of the second kind
% out = dbesselk(l, z)
% Returns d/dz [K_l(z)] using the recurrence formula:
% K_l'(z) = -0.5 * (K_{l-1}(z) + K_{l+1}(z))
out = -0.5 * (besselk(l-1, z) + besselk(l+1, z));
end
function out = dbesselj(l, z)
%DBESSELJ Derivative of the Bessel function of the first kind
% out = dbesselj(l, z)
% Returns d/dz [J_l(z)] using the recurrence formula:
% J_l'(z) = 0.5 * (J_{l-1}(z) – J_{l+1}(z))
out = 0.5 * (besselj(l-1, z) – besselj(l+1, z));
end
function out = dbesseli(l, z)
%DBESSELI Derivative of the modified Bessel function of the first kind
% out = dbesseli(l, z)
% Returns d/dz [I_l(z)] using the recurrence formula:
% I_l'(z) = 0.5 * (I_{l-1}(z) + I_{l+1}(z))
out = 0.5 * (besseli(l-1, z) + besseli(l+1, z));
end
function out = dbesselh(l, z)
%DBESSELH Derivative of the Hankel function of the first kind
% out = dbesselh(l, z)
% Returns d/dz [H_l^{(1)}(z)] using the recurrence formula:
% H_l^{(1)’}(z) = 0.5 * (H_{l-1}^{(1)}(z) – H_{l+1}^{(1)}(z))
out = 0.5 * (besselh(l-1, 1, z) – besselh(l+1, 1, z));
end
function x = bessel0j(l,q,opt)
% a row vector of the first q roots of bessel function Jl(x), integer order.
% if opt = ‘d’, row vector of the first q roots of dJl(x)/dx, integer order.
% if opt is not provided, the default is zeros of Jl(x).
% all roots are positive, except when l=0,
% x=0 is included as a root of dJ0(x)/dx (standard convention).
%
% starting point for for zeros of Jl was borrowed from Cleve Moler,
% but the starting points for both Jl and Jl’ can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13.
%
% David Goodmanson
%
% x = bessel0j(l,q,opt)
k = 1:q;
if nargin==3 && opt==’d’
beta = (k + l/2 – 3/4)*pi;
mu = 4*l^2;
x = beta – (mu+3)./(8*beta) – 4*(7*mu^2+82*mu-9)./(3*(8*beta).^3);
for j=1:8
xnew = x – besseljd(l,x)./ …
(besselj(l,x).*((l^2./x.^2)-1) -besseljd(l,x)./x);
x = xnew;
end
if l==0
x(1) = 0; % correct a small numerical difference from 0
end
else
beta = (k + l/2 – 1/4)*pi;
mu = 4*l^2;
x = beta – (mu-1)./(8*beta) – 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);
for j=1:8
xnew = x – besselj(l,x)./besseljd(l,x);
x = xnew;
end
end
end
% — Local helper function for derivative of Bessel function —
function dJ = besseljd(l, x)
dJ = 0.5 * (besselj(l – 1, x) – besselj(l + 1, x));
end
function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi)
% n = 0 mode
term1 = (d_vec(1) * Z0d * besselj(l, wk*r)) / (dbesselj(l, wk*RI));
% n >= 2 evanescent sum
sum1 = 0;
for nidx = 2:N
sum1 = sum1 + d_vec(nidx) * Znd(nidx) * besseli(l, k(nidx)*r) …
/( dbesseli(l, k(nidx)*RI));
end
% q = 1..Q radial sum
sum2 = 0;
for qidx = 1:Q
sum2 = sum2 + e_vec(qidx) *chi(qidx)* besselj(l, chi(qidx)*r) / dbesselj(l, chi(qidx)*RI);
end
v_D_val = term1 + sum1 + sum2;
end I write the following code but i used the following parametrs and i did not get any numerical instabilities:
RI = 34.5; % Inner radius (m)
ratio = 47.5/34.5; % Ratio R_E / R_I
RE = ratio*RI; % outer radius (m)
h = 200/34.5 * RI; % Water depth (m)
d = 38/34.5 * RI; % Interface depth (m)
b = h-d;
But when i switch the parameters with this then i get numerical instabilities
RI = 0.4; % Inner radius (m)
RE = 0.8; % Outer radius (m)
d = 0.3; % Draft (m)
h = 1.0; % Water depth (m)
b = h-d; %from cylinder bottom to seabed (m)
the following is code
clear all;
close all;
clc;
tic;
% diary(‘iter_log.txt’);
% diary on
% fprintf(‘Run started: %sn’, datestr(now));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters from your setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 10; % Number of N modes
M = 10; % Number of M modes
Q = 10; % Number of Q modes
RI = 0.4; % Inner radius (m)
RE = 0.8; % Outer radius (m)
d = 0.3; % Draft (m)
h = 1.0; % Water depth (m) % Interface depth (m)
b = h-d; % Distance from cylinder bottom to seabed (m)
g = 9.81; % Gravity (m/s^2)
tau =0.2; % porosity ratio
gamma = 1.0; % discharge coefficient
b1 = (1 – tau) / (2 * gamma * tau^2); % nonlinear coeff
tol = 1e-4; % convergence tolerance
max_iter = 100; % max iterations
l = 0; % Azimuthal order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === Compute Bessel roots and scaled chi values ===
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
roots = bessel0j(l, Q); % Zeros of J_l(x)
chi = roots ./ RI; % chi_q^l = roots scaled by radius
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ===== Define Range for k_0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X_kRE = linspace(0.0025, 0.16,100); %%%%%% this is now k_0
X_kRE = linspace(0.05, 4.5, 40);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%===============Initialize omega array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta0 = zeros(size(X_kRE)); % Preallocate
omega = zeros(size(X_kRE)); %%%% omega preallocate
% iters_used = zeros(length(X_kRE),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for idx = 1:length(X_kRE)
wk = X_kRE(idx); % dimensional wavenumber
omega(idx) = sqrt(g * wk * tanh(wk * h)); % dispersion relation
% fprintf(‘n— idx %3d (T=%6.3f s) —n’, idx, 2*pi/omega(idx));
% drawnow;
%%%%%%%%%%%%% Compute group velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cg(idx) = (g * tanh(wk * h) + g * wk * h * (sech(wk * h))^2) …
* omega(idx) / (2 * g * wk * tanh(wk * h));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% =======Compute a1 based on tau
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a1(idx) = 0.93 * (1 – tau) * Cg(idx);
% dissip(idx) = a1(idx)/omega(idx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%=============== Find derivative of Z0 at z = -d
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z0d = (cosh(wk * h) * wk * sinh(wk * b)) / (2 * wk * h + sinh(2 * wk * h));
fun_alpha = @(x) omega(idx)^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes
opts_lsq = optimoptions(‘lsqnonlin’,’Display’,’off’); % ← add this
for n = 1:N
if n == 1
k(n) = -1i * wk;
else
x0_left = (2*n – 3) * pi / (2*h);
x0_right = (2*n – 1) * pi / (2*h);
guess = (x0_left + x0_right)/2;
% Use lsqnonlin for better convergence
k(n) = lsqnonlin(fun_alpha, guess, x0_left, x0_right,opts_lsq);
%%%%%%%%%%%%% Derivative of Z_n at z = -d %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Znd(n) = -k(n) * (cos(k(n)*h) * sin(k(n)*b)) / (2*k(n)*h + sin(2*k(n)*h));
end
end
%% finding the root lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;
end
% %%%%%%%%%%%%%%%%%%%%% Define matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = zeros(M,N);
A(1,1) = (cosh(wk*h)*sinh(wk*b))/(wk*b*(2*wk*h+sinh(2*wk*h)))*(besselh(l, wk*RE)) /(dbesselh(l, wk*RE)); % Set A_00
for j = 2:N
A(1,j) =(cos(k(j)*h)*sin(k(j)*b) )/ (k(j)*b*(2*k(j)*h + sin(2*k(j)*h)))*( besselk(l, k(j)*RE) )…
/(dbesselk(l, k(j)*RE)); % Set A_0n
end
for i = 2:M
A(i,1) =(2 * cosh(wk*h) * (-1)^(i-1) * wk*sinh(wk*b)) / (b * (2*wk*h + sinh(2*wk*h))*(wk^2 + m(i)^2))…
* (besselh(l, wk*RE))/ (dbesselh(l, wk*RE)); % Set A_m0
end
for i = 2:M
for j = 2:N
A(i,j) = (2*cos(k(j)*h)*(-1)^(i-1) *k(j)* sin(k(j)*b))/ (b * (2*k(j)*h + sin(2*k(j)*h))*(k(j)^2-m(i)^2))…
*(besselk(l, k(j)*RE)) / (dbesselk(l, k(j)*RE));% Set A_mn
end
end
%%%%%%%%%%%%%%%%%%%%%%%% DefineMatrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N,M);
B(1,1) = (4*sinh(wk*b)) / (RE * wk*log(RE/RI) *cosh(wk*h)); %set B_00
%%%%%%%%%%%B_0m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(1,j) = Rml_prime_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%%%%%%%%B_n0%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
B(i,1)=(4*sin(k(i)*b))/(RE*k(i)*log(RE/RI)*cos(k(i)*h));% Set B_n0
end
%%%%%%%%%%%%%%%%%%%%%%%%%%B_nm%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
for j=2:M
Rml_prime_RE = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
B(i,j) = Rml_prime_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * wk*log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for j = 2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(1,j) = Rml_prime_star_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
% %%%%%%% Define Cn0%%%%%%
for i = 2:N
C(i,1) = -4 * sin(k(i)*b) / (RE *k(i)* log(RE/RI) * cos(k(i)*h));
end
% % %%%%%% Define Cnm%%%%%%
for i = 2:N
for j =2:M
Rml_prime_star_RE = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))…
/(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
C(i,j) = Rml_prime_star_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%% write Matrix D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = zeros(M,N);
%%% write D00 %%%%%%
D(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2*wk*h + sinh(2*wk*h))) * (besselj(l, wk*RI))/ (dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for j =2:N
D(1,j) = (cos(k(j)*h) * sin(k(j)*b)) / (k(j)*b * (2*k(j)*h + sin(2*k(j)*h))) * (besseli(l, k(j)*RI) )/ (dbesseli(l, k(j)*RI));
end
%%%%%% write Dm0 %%%%%%%%%
for i = 2:M
D(i,1) = (2 * cosh(wk*h) * (-1)^(i-1) * wk * sinh(wk*b)) /(b * (2*wk*h + sinh(2*wk*h)) * (wk^2 + m(i)^2))…
*(besselj(l, wk*RI) )/(dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for i = 2:M
for j = 2:N
D(i,j) = (2 * cos(k(j)*h) * (-1)^(i-1) * k(j) * sin(k(j)*b)) /(b * (2*k(j)*h + sin(2*k(j)*h)) * (k(j)^2 – m(i)^2))…
*(besseli(l, k(j)*RI)) / (dbesseli(l, k(j)*RI));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M); % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
E(1,1) = (4 * sinh(wk*b)) / (RI *wk* log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(1,j) = Rml_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
%
end
%%%%%% Define Eno%%%%%%%%%%%%%%
for i = 2:N
E(i,1) = (4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
for i = 2:N
for j = 2:M
Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) …
– besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))…
/ (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(i,j) = Rml_prime_RI * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * wk*log(RE/RI) * cosh(wk*h));
% %%%%%% Define F0m%%%%
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(1,j) = Rml_star_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%% Defin Fn0%%%%%%
for i = 2:N
F(i,1) = (-4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
%%%%%%% Define Fnm %%%%%%%
for i = 2:N
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) …
– besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) …
* besseli(l, m(j)*RE) – besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(i,j) = Rml_star_prime_RI * (4 * k(i)* (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 – m(j)^2));
end
end
% % %%%%%%%%%% Write Matrix W %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = zeros(N,Q);
% %%W_{0q}
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 – r); % coth(chi*b) via decaying exp
W(1,q) = -(4*chi(q)/cosh(wk*h))*(( wk*sinh(wk*b)*coth_chib – chi(q)*cosh(wk*b) ) / ( wk^2 – chi(q)^2 ) ); %prof.chen
%
% % W(1, q) = – (4 * chi(q) / cosh(wk * h)) * (chi(q) * sinh(chi(q) * b) * cosh(wk * b) – wk * sinh(wk * b) * cosh(chi(q) * b)) …
% % / ((chi(q)^2 – wk^2) * sinh(chi(q) * b));
end
% % %%% Write W_{nq}
for n = 2:N
for q = 1:Q
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
coth_chib = (1 + r) / (1 – r); % coth(chi*b), stable
W(n,q) = -(4*chi(q)/cos(k(n)*h)) * …
( ( k(n)*sin(k(n)*b)*coth_chib + chi(q)*cos(k(n)*b) ) … %prof.chen
/ ( k(n)^2 + chi(q)^2 ) );
% % W(n, q) = – (4 * chi(q) / cos(k(n) * h)) * (chi(q) * sinh(chi(q) * b) * cos(k(n) * b) + k(n) * sin(k(n) * b) * cosh(chi(q) * b)) …
% % / ((chi(q)^2 + k(n)^2) * sinh(chi(q) * b));
end
end
%%%%%%%%%%%%%%%%%%%%%%%Find E1_{q}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f2 = omega(idx)^2 / g;
E1 = zeros(Q,1);
S_q = zeros(Q,1);
for q = 1:Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% this expression is same as prof. Chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
exp2chi_d = exp(-2 * chi(q) * d); % e^{-2*chi*d}
exp2chi_b = exp(-2 * chi(q) * b); % e^{-2*chi*(h-d)}
exp2chi_h = exp(-2 * chi(q) * h); % e^{-2*chi*h}
exp1chi_d = exp(-chi(q) * d); % e^{-chi*d}
% intermediate ratio (same order as Fortran)
BB = (f2 – chi(q)) * (exp2chi_b / exp2chi_d);
BB = BB + (f2 + chi(q)) * (1 + 2 * exp2chi_b);
BB = BB / (f2 – chi(q) + (f2 + chi(q)) * exp2chi_d);
BB = BB / (1 + exp2chi_b);
% final coefficient
E1(q) = 2 * (BB * exp2chi_d – 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This expression is that one derived directly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% % N0, D0, H0 exactly as in the math above
% N0 = (chi(q) – f2) + (chi(q) + f2) * exp2chi_h;
% D0 = (f2 – chi(q)) + (f2 + chi(q)) * exp2chi_d;
% H0 = 1 + exp2chi_b;
%
% % final coefficient (no add/subtract step)
% E1(q) = 2 * N0 / (D0 * H0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This one is the hyperbolic one form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% num = chi(q)*cosh(chi(q)*h) – f2*sinh(chi(q)*h);
% den = (f2*cosh(chi(q)*d) – chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*(h – d)); % or: * cosh(chi(q)*b)
% E1(q) = num / den;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% This is the upper region velocity potential at z=0, write
%%%%%%%%%%%%%%%% from prof.chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CCmc = 1 + (exp2chi_b / exp2chi_d) + 2*exp2chi_b;
CCmc = CCmc / (1 + exp2chi_b);
CCmc = BB * (1 + exp2chi_d) – CCmc;
S_q(q) = CCmc * exp1chi_d; % this equals the summand value for this q
% — direct hyperbolic (exact) —
% % % C_hyp = ( chi(q)*cosh(chi(q)*h) – f2*sinh(chi(q)*h) ) …
% % % / ( (f2*cosh(chi(q)*d) – chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*b) );
% % %
% % % S_hyp(q) = C_hyp * cosh(chi(q)*d) + sinh(chi(q)*h) / cosh(chi(q)*b); % summand
end
% H(q): diagonal term (exponential form)
H1 = zeros(Q,1);
for q = 1:Q
% dissip =
r = exp(-2*chi(q)*b); % e^{-2*chi*b}
mid = 4*r/(1 – r^2) – E1(q); % 2/sinh(2*chi*b) – C_m
H1(q) = mid + 1i*a1(idx)*chi(q)/omega(idx); % mid + i*(a1/w)*chi
end
H = diag(H1);
%%%%%%%%%%matric G%%%%%%%%%%%%%%%%%%
G = zeros(Q, N); % Preallocate
for q = 1:Q
G(q,1) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Z0d …
* ( chi(q) / (chi(q)^2 – wk^2) )* ( besselj(l, wk*RI) / dbesselj(l, wk*RI) ); % this one is my calculated
end
% %G_qn
for q = 1:Q
for n = 2:N
G(q, n) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Znd(n) …
* ( chi(q) / (k(n)^2 + chi(q)^2) )* ( besseli(l, k(n)*RI) / dbesseli(l, k(n)*RI) ); % this one is my calculated
end
end
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Define the right hand side vector
U = zeros(2*M + 2*N + Q, 1); % Full vector, M+N+M+N+Q elements
% % Block 1: U (size M = 4)
for i = 1:M
if i == 1
U(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
U(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2)*cosh(wk*h))*besselj(l, wk*RE); %Z_m^l
end
end
% % Block 2: Y (size N = 3)
for j = 1:N
if j == 1
U(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
U(j + M, 1) = 0; % Y_n^l
end
end
% % Block 3: X (size M = 4)
for i = 1:M
U(i + M + N, 1) = 0; % X_0^l, X_m^l
end
% Block 4: W (size N = 3)
for j = 1:N
U(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
% % Identity matrices
I_M = eye(M);
I_N = eye(N);
% Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);
ZMQ = zeros(M,Q);
ZNQ = zeros(N,Q);
ZQM = zeros(Q,M);
ZQN = zeros(Q,N);
% Construct the full block matrix S
S = [ I_M, -A, ZMM, ZMN, ZMQ;
-B, I_N, -C, ZNN, ZNQ;
ZMM, ZMN, I_M, -D, ZMQ;
-E, ZNN, -F, I_N, -W;
ZQM, ZQN, ZQM, -G, H ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
converged = false;
psi = zeros(Q,1);
T_old = []; % <– track the full solution from prev. iter
% T_old = zeros(2*M + 2*N + Q, 1); % previous unknowns (seed)
for iter = 1:max_iter
% (A) overwrite ONLY Block 5 of U from current psi
for q = 1:Q
U(2*M + 2*N + q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi(q);
end
T = S U; % Solve the linear system
% T = pinv(S) * U; % Use pseudoinverse for stability
b_vec = T(1:M); % Coefficients b^l
a_vec = T(M+1 : M+N); % Coefficients a^l
c_vec = T(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = T(2*M+N+1:2*M+2*N); % Coefficients d^l
e_vec = T(2*M+2*N+1:end); % (Q×1)
% (D) update psi for NEXT iteration from CURRENT coefficients
for q = 1:Q
integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) …
.* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) …
.* besselj(l, chi(q)*r) .* r;
psi(q) = integral(integrand, 0, RI, ‘AbsTol’,1e-8, ‘RelTol’,1e-6);
end
% % % drawnow; % optional: flush output each iter
% === compact per-iteration print, only for T in [5,35] ===
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf(‘idx %3d | iter %2d | ||psi|| = %.3en’, idx, iter, norm(psi));
% drawnow;
% end
% % % % (E) convergence on ALL unknowns (T vs previous T)
% — per-iteration diagnostics (optional) —
if ~isempty(T_old)
diff_T = max(abs(T – T_old));
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
% fprintf(‘k0 idx %3d | iter %2d: max|ΔT| = %.3e, ||psi|| = %.3en’, …
% idx, iter, diff_T, norm(psi));
% drawnow; % ← add this to flush each iteration line too
% end
if diff_T < tol
converged = true;
break
end
end
T_old = T;
end % <<< end of iter loop
% % —- ONE summary line per frequency (place it here) —-
% iters_used(idx) = iter; % how many iterations this frequency needed
% fprintf(‘idx %3d (T=%6.3f s): iters=%2dn’, idx, 2*pi/omega(idx), iter);
% drawnow;
% ————– end nonlinear iteration ————–
% % %%%%%%%%% Wave motion%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
term1 = (d_vec(1)*cosh(wk*h)^2) / ((2*wk*h + sinh(2*wk*h)) * dbesselj(l, wk*RI));
sum1 = 0;
for n =2:N
sum1 = sum1 + d_vec(n)*cos(k(n)*h)^2/ ((2*k(n)*h + sin(2*k(n)*h))* dbesseli(l, k(n)*RI));
end
phiUell = 0;
for q =1:Q
phiUell = phiUell +e_vec(q)* S_q(q);
end
eta0(idx) = abs(term1+sum1+phiUell);
end
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % Plotting Data
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% figure(2); clf;
% >>> CHANGES START (3): one figure = two taus + Dr. Chen only
% plot(2*pi./omega, eta0, ‘k’, ‘LineWidth’, 1.5); %%% this is T which is T= 2*pi/omega plotting against T
plot(omega.^2 * RE / g, eta0, ‘k’, ‘LineWidth’, 1.5); %%% this is T which is T= 2*pi/omega plotting against T
hold on; % Add extra plots on same figure
% % % %
% % % % % Now plot the 3 CSV experimental points
% scatter(data_chen(:,1), data_chen(:,2), 30, ‘r’, ‘o’, ‘filled’);
% % scatter(data_liu(:,1), data_liu(:,2), 30, ‘g’, ‘s’, ‘filled’);
% scatter(data_exp(:,1), data_exp(:,2), 30, ‘b’, ‘^’, ‘filled’);
xlabel(‘$T$’, ‘Interpreter’, ‘latex’);
ylabel(‘$|eta / (iA)|$’, ‘Interpreter’, ‘latex’);
title(‘Wave motion amplitude for $R_E = 138$’, ‘Interpreter’, ‘latex’);
legend({‘$tau=0.2$’,’Model test’}, …
‘Interpreter’,’latex’,’Location’,’northwest’);
grid on;
% xlim([5 35]); % Match the reference plot
% ylim([0 4.5]); % Optional: based on expected peak height
xlim([0 4]); % Match the reference plot
ylim([0 7]); % Optional: based on expected peak height
% === plotting section ===
% diary off
elapsedTime = toc;
disp([‘Time consuming = ‘, num2str(elapsedTime), ‘ s’]);
function out = dbesselk(l, z)
%DBESSELK Derivative of the modified Bessel function of the second kind
% out = dbesselk(l, z)
% Returns d/dz [K_l(z)] using the recurrence formula:
% K_l'(z) = -0.5 * (K_{l-1}(z) + K_{l+1}(z))
out = -0.5 * (besselk(l-1, z) + besselk(l+1, z));
end
function out = dbesselj(l, z)
%DBESSELJ Derivative of the Bessel function of the first kind
% out = dbesselj(l, z)
% Returns d/dz [J_l(z)] using the recurrence formula:
% J_l'(z) = 0.5 * (J_{l-1}(z) – J_{l+1}(z))
out = 0.5 * (besselj(l-1, z) – besselj(l+1, z));
end
function out = dbesseli(l, z)
%DBESSELI Derivative of the modified Bessel function of the first kind
% out = dbesseli(l, z)
% Returns d/dz [I_l(z)] using the recurrence formula:
% I_l'(z) = 0.5 * (I_{l-1}(z) + I_{l+1}(z))
out = 0.5 * (besseli(l-1, z) + besseli(l+1, z));
end
function out = dbesselh(l, z)
%DBESSELH Derivative of the Hankel function of the first kind
% out = dbesselh(l, z)
% Returns d/dz [H_l^{(1)}(z)] using the recurrence formula:
% H_l^{(1)’}(z) = 0.5 * (H_{l-1}^{(1)}(z) – H_{l+1}^{(1)}(z))
out = 0.5 * (besselh(l-1, 1, z) – besselh(l+1, 1, z));
end
function x = bessel0j(l,q,opt)
% a row vector of the first q roots of bessel function Jl(x), integer order.
% if opt = ‘d’, row vector of the first q roots of dJl(x)/dx, integer order.
% if opt is not provided, the default is zeros of Jl(x).
% all roots are positive, except when l=0,
% x=0 is included as a root of dJ0(x)/dx (standard convention).
%
% starting point for for zeros of Jl was borrowed from Cleve Moler,
% but the starting points for both Jl and Jl’ can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13.
%
% David Goodmanson
%
% x = bessel0j(l,q,opt)
k = 1:q;
if nargin==3 && opt==’d’
beta = (k + l/2 – 3/4)*pi;
mu = 4*l^2;
x = beta – (mu+3)./(8*beta) – 4*(7*mu^2+82*mu-9)./(3*(8*beta).^3);
for j=1:8
xnew = x – besseljd(l,x)./ …
(besselj(l,x).*((l^2./x.^2)-1) -besseljd(l,x)./x);
x = xnew;
end
if l==0
x(1) = 0; % correct a small numerical difference from 0
end
else
beta = (k + l/2 – 1/4)*pi;
mu = 4*l^2;
x = beta – (mu-1)./(8*beta) – 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);
for j=1:8
xnew = x – besselj(l,x)./besseljd(l,x);
x = xnew;
end
end
end
% — Local helper function for derivative of Bessel function —
function dJ = besseljd(l, x)
dJ = 0.5 * (besselj(l – 1, x) – besselj(l + 1, x));
end
function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi)
% n = 0 mode
term1 = (d_vec(1) * Z0d * besselj(l, wk*r)) / (dbesselj(l, wk*RI));
% n >= 2 evanescent sum
sum1 = 0;
for nidx = 2:N
sum1 = sum1 + d_vec(nidx) * Znd(nidx) * besseli(l, k(nidx)*r) …
/( dbesseli(l, k(nidx)*RI));
end
% q = 1..Q radial sum
sum2 = 0;
for qidx = 1:Q
sum2 = sum2 + e_vec(qidx) *chi(qidx)* besselj(l, chi(qidx)*r) / dbesselj(l, chi(qidx)*RI);
end
v_D_val = term1 + sum1 + sum2;
end bessel functions, system of equations MATLAB Answers — New Questions
Checking Where Tenant Users Go as Guests
Managing External Guest Activity Controlled by the Host Tenant
After the original fuss about the Teams feature where users can invite external people to join chats through email had subsided, a brief flame reignited when The Hacker News posted an article explaining the risks involved when someone invites a user from your tenant to join them for a chat. If the invitation is accepted, a guest account for the invited user is created in the host tenant and the chat proceeds.
The problem highlighted by the article is that the host tenant determines what kind of protection it wishes to apply to Teams messaging. The host tenant can turn Microsoft Defender for Office 365 protection off, and disable the protections against weaponized file types and malicious URLs recently introduced for Teams. With all protection disabled, there’s nothing to stop a malicious file or URL being shared in chat that could potentially infect your tenant.
It’s correct that Entra B2B Collaboration works this way. When someone becomes a guest in another tenant, they agree to collaboration managed by the host tenant rather than their home tenant. For instance, any compliance records captured to note actions taken by the guest account are owned by the host tenant. The user’s home tenant has no idea whatsoever about what someone does outside the security boundary of that tenant.
Discovering Where Tenant Users Go as Guests
This then begs the question whether a home tenant can discover what host tenants their users are active in. A home tenant can restrict access to certain host tenants through policy, but that doesn’t mean that guest accounts for users in the home tenant don’t exist and are not active in “blocked” host tenants. If those accounts were created prior to the block being imposed, they will continue to function afterwards.
Microsoft’s solution for figuring out where users are active guests is a workbook to report cross-tenant activity, which uses sign-in data ingested into Azure Log Analytics to generate an activity analysis. This is great if the tenant ingests Entra ID sign-in logs into Log Analytics because you can go back further than the 30-day retention period for online sign-in logs used by Entra ID. Unfortunately, not every tenant uses Log Analytics, but we can use PowerShell to do the job. Let’s look at what needs to be done.
Using PowerShell to Analyze External Guest Activity
PowerShell can access the same sign-in audit data using cmdlets from the Microsoft Graph PowerShell SDK. In this case, we use the Get-MgBetaAuditLogSignIn cmdlet to fetch sign-in logs because the beta version delivers more information than the production (V1.0) version does. Two key properties are the HomeTenantId and the ResourceTenantId. The first is the identifier for the home tenant; the second is the identifier for the host tenant. For normal activity by a user within their home tenant, the two values are the same. When the identifiers differ, we know we’re dealing with cross-tenant activity.
The logic flow is straightforward:
- Find audit log records for successful sign-in records for the period to analyze. A very large number of sign-in records might be available for large tenants, so it will take time to retrieve the audit data for a complete month. Unfortunately, the Graph API doesn’t support a server-side filter to find records where the home tenant identifier is different to the host tenant identifier, so a client-side filter is necessary to extract the records for cross-tenant activity.
- Cross-tenant activity can be outbound (our users accessing host tenants as guests) or inbound (users from other tenants accessing our tenant as guests). The script creates separate arrays for both types of activity.
- For each tenant identifier, the script calls the Find-MgTenantRelationshipTenantInformationByTenantId cmdlet to resolve the identifier to a tenant name and domain. For example, the identifier 72f988bf-86f1-41af-91ab-2d7cd011db47 resolves to “Microsoft.”
- The script extracts details of the number of unique sign-ins for the tenant, the user principal names of the accounts that are signing in as guests, and apps hosted by the tenant that are being signed into.
- The script reports some basic information and generates an output file. If the ImportExcel module is available, the output is an Excel worksheet. Otherwise, it is a CSV file.
Figure 1 shows the external tenant access report for my tenant. I’m the only one who has accessed other tenants as a guest, but a reasonable number of people from other tenants have guest accounts in my tenant, notably to help with the Office 365 for IT Pros eBook.
The script is available from the Office 365 for IT Pros GitHub repository.
Only Sign-ins
It’s important to recognize that the data reported here are sign-in records captured by Entra ID when someone goes through the sigh-in process. There’s no detail about what someone actually does after they successfully authenticate and connect to a host tenant. That data resides in the host tenant and is inaccessible unless an administrator of that tenant makes it available to you. But at least we can determine who’s an active guest in other tenants, which is the purpose of this exercise.
Need help to write and manage PowerShell scripts for Microsoft 365, including Azure Automation runbooks? Get a copy of the Automating Microsoft 365 with PowerShell eBook, available standalone or as part of the Office 365 for IT Pros eBook bundle.
Linux ubuntu 22.04 lts videoreader h265
Hello, how can i manage to read a video H265 in matlab on ubuntu 22.04 lts ? If you have an idea on how can I proceed I will be pleased to have your help.Hello, how can i manage to read a video H265 in matlab on ubuntu 22.04 lts ? If you have an idea on how can I proceed I will be pleased to have your help. Hello, how can i manage to read a video H265 in matlab on ubuntu 22.04 lts ? If you have an idea on how can I proceed I will be pleased to have your help. linux, h265, video MATLAB Answers — New Questions
the dimension of matrix when using ‘sparse’ to generate sparse matrix
Why matlab has removed the high dimension operation of sparse?
I used to generate 3D sparse matrix using sparse. However, I find an error message when using this function says: Error using sparse, inputs must be 2-D.
Is this problem solvable?Why matlab has removed the high dimension operation of sparse?
I used to generate 3D sparse matrix using sparse. However, I find an error message when using this function says: Error using sparse, inputs must be 2-D.
Is this problem solvable? Why matlab has removed the high dimension operation of sparse?
I used to generate 3D sparse matrix using sparse. However, I find an error message when using this function says: Error using sparse, inputs must be 2-D.
Is this problem solvable? sparse matrix generation, 3d matrix MATLAB Answers — New Questions
MATLAB Academy Won’t Load
Hi there!! I am trying to get through some of the MATLAB academy online courses, but they won’t open. The tab will open but then it just spins for hours and hours and won’t actually let me work on them. This has been happening intermittently for about a week now and I never had a problem before this. I am using an HP Ominibook 7. I already tried clearing my cache and cookies and restarting my computer, but still the issue persists. Does anybody have any ideas for what’s going on and how to fix it?Hi there!! I am trying to get through some of the MATLAB academy online courses, but they won’t open. The tab will open but then it just spins for hours and hours and won’t actually let me work on them. This has been happening intermittently for about a week now and I never had a problem before this. I am using an HP Ominibook 7. I already tried clearing my cache and cookies and restarting my computer, but still the issue persists. Does anybody have any ideas for what’s going on and how to fix it? Hi there!! I am trying to get through some of the MATLAB academy online courses, but they won’t open. The tab will open but then it just spins for hours and hours and won’t actually let me work on them. This has been happening intermittently for about a week now and I never had a problem before this. I am using an HP Ominibook 7. I already tried clearing my cache and cookies and restarting my computer, but still the issue persists. Does anybody have any ideas for what’s going on and how to fix it? matlab academy, loading course hangs MATLAB Answers — New Questions
Steering System Rack and Pinion
Dear Sir,
I am modeling an steering system so that i can control the steering wheel with DC motor but I am now stuck at modeling steering system. I need the model to acting with steering force. When i use other Steering Block it is just right time respond (i thought that didnot contain steering force in that block).
I am using Steering System, Fiala Wheel 2DOF, Vehicle Rigid Body 3DOF in Vehicle Dynamic Blockset.
Fiala Wheel for 2 front wheel with output are Fx, Fy, Fz, Mx, My, Mz. Iam using this output là 1×12 vector for Tire Feedback Steering System. Then using Torque input for Steering System to apply for Rigid Body 3DOF.
There are some problems:
When i apply 0 AxlTrq to Fiala Wheel 2DOF and 0 Torque input or 0 Angle Input for Steering System, the model is always turning in 1 direction.
Here is what i get from Graph XY.
I dont know how to fix this.
2. Is AxlTrq needed in Fiala Wheel 2Dof? When i apply AxlTrq = Fx*Re, Fx is the output from front right and left force from Vehicle Rigid Body, Re is feedback from Fiala Wheel Radius Effective.
I got a problem:
But I think this can be fixxed if the problem 1 can be solved.
3. I wonder if my model is doing right?
Thank you for your help.Dear Sir,
I am modeling an steering system so that i can control the steering wheel with DC motor but I am now stuck at modeling steering system. I need the model to acting with steering force. When i use other Steering Block it is just right time respond (i thought that didnot contain steering force in that block).
I am using Steering System, Fiala Wheel 2DOF, Vehicle Rigid Body 3DOF in Vehicle Dynamic Blockset.
Fiala Wheel for 2 front wheel with output are Fx, Fy, Fz, Mx, My, Mz. Iam using this output là 1×12 vector for Tire Feedback Steering System. Then using Torque input for Steering System to apply for Rigid Body 3DOF.
There are some problems:
When i apply 0 AxlTrq to Fiala Wheel 2DOF and 0 Torque input or 0 Angle Input for Steering System, the model is always turning in 1 direction.
Here is what i get from Graph XY.
I dont know how to fix this.
2. Is AxlTrq needed in Fiala Wheel 2Dof? When i apply AxlTrq = Fx*Re, Fx is the output from front right and left force from Vehicle Rigid Body, Re is feedback from Fiala Wheel Radius Effective.
I got a problem:
But I think this can be fixxed if the problem 1 can be solved.
3. I wonder if my model is doing right?
Thank you for your help. Dear Sir,
I am modeling an steering system so that i can control the steering wheel with DC motor but I am now stuck at modeling steering system. I need the model to acting with steering force. When i use other Steering Block it is just right time respond (i thought that didnot contain steering force in that block).
I am using Steering System, Fiala Wheel 2DOF, Vehicle Rigid Body 3DOF in Vehicle Dynamic Blockset.
Fiala Wheel for 2 front wheel with output are Fx, Fy, Fz, Mx, My, Mz. Iam using this output là 1×12 vector for Tire Feedback Steering System. Then using Torque input for Steering System to apply for Rigid Body 3DOF.
There are some problems:
When i apply 0 AxlTrq to Fiala Wheel 2DOF and 0 Torque input or 0 Angle Input for Steering System, the model is always turning in 1 direction.
Here is what i get from Graph XY.
I dont know how to fix this.
2. Is AxlTrq needed in Fiala Wheel 2Dof? When i apply AxlTrq = Fx*Re, Fx is the output from front right and left force from Vehicle Rigid Body, Re is feedback from Fiala Wheel Radius Effective.
I got a problem:
But I think this can be fixxed if the problem 1 can be solved.
3. I wonder if my model is doing right?
Thank you for your help. vehicle dynamic blockset, steering system, rack and pinion, simulink, fiala wheel 2dof, vehicle rigid body 3dof MATLAB Answers — New Questions
Microsoft Increases Office 365 and Microsoft 365 License Prices
New Microsoft 365 Pricing Goes into Effect on July 1, 2026
On December 4, 2025, Microsoft announced a range of price increases for Microsoft 365 monthly licenses. The new pricing (Figure 1) goes into effect from July 1, 2026, the start of Microsoft’s FY27 fiscal year.
According to Microsoft, they want to “give customers ample time to plan.” However, there’s not much choice for tenants if their operations are embedded in the Microsoft 365 ecosystem, so this is a case of “getting used to new pricing” rather than “having time to consider migrating away from Microsoft 365.” Once you’re embedded in the Microsoft 365 ecosystem, it’s hard to leave.
Some organizations do consider going back to on-premises servers. It’s certainly an option, even to the now available and oddly named Microsoft 365 Local, a product that shares precisely nothing but its name with the rest of the Microsoft 365 ecosystem.
Last Microsoft 365 License Increase in 2022
Microsoft last increased Microsoft 365 license prices in March 2022. At the time, Microsoft added $3/monthly to Office 365 E3m and E5, and $4/monthly to Microsoft 365 E3. The Microsoft 365 E5 price was left unchanged.
This time round, the monthly increases range from zero (Office 365 E1) to $3 (the big plans used by large enterprises like Office 365 E3 and Microsoft 365 E5). At $2/average across the Microsoft 365 base (around 446 million paid seats based on data provided at Microsoft’s FY26 Q1 earnings), the increase could bring in an extra $10.7 billion. The price changes shown in Figure 1 apply to the commercial cloud. Equivalent increases apply to other sectors, such as education and government.
In FY26 Q1, the Microsoft Cloud operated at a healthy 68% operating margin, so it’s not as if Microsoft does not achieve an adequate return from Microsoft 365. However, as noted in the earnings transcript, the operating margin for the Microsoft Cloud is down year-over-year due to “investments in AI.” One interpretation is that the extra $10 billion from the price increases will offset some of the red ink Microsoft is bleeding because of the investments they’re making in datacenter capacity, hardware, and software needed to make Copilot useful,
Justifying the Additional Cost
Just like last time around, Microsoft justifies the increase by pointing to an array of new features and functionality that they’ve delivered. Microsoft 365 E5 customers recently received news that they will soon get Security Copilot, and another announcement revealed that the Microsoft 365 E3 and E5 plans will both gain functionality from the Microsoft Intune Suite in the coming months.
Plans that don’t get Security Copilot or the Intune Suite must do with new apps like Microsoft Loop, Clipchamp, and Places, all introduced since the 2022 price increase. Good as these apps are, a tenant has to use them to extract value to justify the additional cost,. A welcome change is the addition of Microsoft 365 Defender for Office 365 P1 to Office 365 E3 and Microsoft 365 E3, even if this might provoke further worry about incurring cost to license shared mailboxes that benefit from Defender functionality.
So Many New Features
Curiously, the blog highlights the release of 1,100 new features in the last year across “Microsoft 365, Copilot, and SharePoint.” I thought SharePoint was a core part of Microsoft 365, but apparently, it’s so important that SharePoint deserves its own mention. Teams just doesn’t get a mention these days. I also wonder how many of the new features are related to Copilot and are therefore useless to tenants that don’t use Copilot.
By comparison, in 2022, Microsoft claimed the release of 1,400 new features in communication and collaboration (aka Teams), security and compliance, and AI and automation (not Copilot!). At the time, I asked how many of the updates were useful. The same could be asked now. Quantity of updates pushed out in a never-ending stream is no substitute for usefulness or quality.
A Question of Value
I’m unsure if any organization can use all the functionality bundled into Microsoft 365. It’s a feature-rich environment with lots to recommend it. I worry about quality of software, the pace of change, the way that Microsoft relentlessly pushes AI at every opportunity, and poor communication about the value of changes at times.
Overall, Microsoft 365 remains very a competitive offering, even if the basic enterprise license is now $312/user/year and the headline E5 license a whopping $720/user/year. Then again, it wasn’t too long ago since a shrink-wrapped copy of Office cost over $300, so perhaps the cost isn’t so bad after all. Either way, I’m sure the increases will cause tenants to devote some time to study their current license mix and allocation to see if any savings are possible (the Microsoft 365 licensing report script might be useful here).
Support the work of the Office 365 for IT Pros team by subscribing to the Office 365 for IT Pros eBook. Your support pays for the time we need to track, analyze, and document the changing world of Microsoft 365 and Office 365. Only humans contribute to our work!
MathWorksSimulation Plugin for Unreal Engine could not be loaded.
Hello,
I am having trouble using the MathWorksSimulation plugin in Unreal Engine for Unreal-Simulink cosimulation. I followed the steps detailed at the link below.
https://www.mathworks.com/help/sl3d/customize-scenes-in-unreal-engine-for-3d-simulations.html
When I open a project an Unreal Engine 5.3 and attempt to add the plugin, it prompts me to reload the scene. I hit yes to rebuild, and Unreal opens the scene is Visual Studio 2022, which is properly configured according to the Unreal Recommendations at the link below.
https://dev.epicgames.com/documentation/en-us/unreal-engine/setting-up-visual-studio-development-environment-for-cplusplus-projects-in-unreal-engine?application_version=5.3
Once I have rebuilt in Visual Studio 2022, and attempt to re-launch the project, I get this error:
"Plugin ‘MathWorksSimulation’ failed to load because module ‘MathWorksSimulation’ could not be loaded. There may be an operating system error or the module may not be properly set up."
Does anyone know how to resolve this? My computer meets the hardware and software specifications detailed here. All the files are in the correct location for Unreal to access them as far as I can verify.
Thanks in advance!
KoryHello,
I am having trouble using the MathWorksSimulation plugin in Unreal Engine for Unreal-Simulink cosimulation. I followed the steps detailed at the link below.
https://www.mathworks.com/help/sl3d/customize-scenes-in-unreal-engine-for-3d-simulations.html
When I open a project an Unreal Engine 5.3 and attempt to add the plugin, it prompts me to reload the scene. I hit yes to rebuild, and Unreal opens the scene is Visual Studio 2022, which is properly configured according to the Unreal Recommendations at the link below.
https://dev.epicgames.com/documentation/en-us/unreal-engine/setting-up-visual-studio-development-environment-for-cplusplus-projects-in-unreal-engine?application_version=5.3
Once I have rebuilt in Visual Studio 2022, and attempt to re-launch the project, I get this error:
"Plugin ‘MathWorksSimulation’ failed to load because module ‘MathWorksSimulation’ could not be loaded. There may be an operating system error or the module may not be properly set up."
Does anyone know how to resolve this? My computer meets the hardware and software specifications detailed here. All the files are in the correct location for Unreal to access them as far as I can verify.
Thanks in advance!
Kory Hello,
I am having trouble using the MathWorksSimulation plugin in Unreal Engine for Unreal-Simulink cosimulation. I followed the steps detailed at the link below.
https://www.mathworks.com/help/sl3d/customize-scenes-in-unreal-engine-for-3d-simulations.html
When I open a project an Unreal Engine 5.3 and attempt to add the plugin, it prompts me to reload the scene. I hit yes to rebuild, and Unreal opens the scene is Visual Studio 2022, which is properly configured according to the Unreal Recommendations at the link below.
https://dev.epicgames.com/documentation/en-us/unreal-engine/setting-up-visual-studio-development-environment-for-cplusplus-projects-in-unreal-engine?application_version=5.3
Once I have rebuilt in Visual Studio 2022, and attempt to re-launch the project, I get this error:
"Plugin ‘MathWorksSimulation’ failed to load because module ‘MathWorksSimulation’ could not be loaded. There may be an operating system error or the module may not be properly set up."
Does anyone know how to resolve this? My computer meets the hardware and software specifications detailed here. All the files are in the correct location for Unreal to access them as far as I can verify.
Thanks in advance!
Kory unreal engine, unreal engine plugin, mathworkssimulation MATLAB Answers — New Questions
0Hz Peak and Absurd Results in Spectrum Analyzer
I wanted to measure the THD of my circuit in MATLAB. I ran a few simulations but got absurd results. To understand, I built a test circuit to see what was happening. It observed a 0 Hz peak and an absurdly high level of THD. I’m not using any code just simscape models, so I don’t understand why this occurs. How can I fix this issue? I also attached the screenshots.I wanted to measure the THD of my circuit in MATLAB. I ran a few simulations but got absurd results. To understand, I built a test circuit to see what was happening. It observed a 0 Hz peak and an absurdly high level of THD. I’m not using any code just simscape models, so I don’t understand why this occurs. How can I fix this issue? I also attached the screenshots. I wanted to measure the THD of my circuit in MATLAB. I ran a few simulations but got absurd results. To understand, I built a test circuit to see what was happening. It observed a 0 Hz peak and an absurdly high level of THD. I’m not using any code just simscape models, so I don’t understand why this occurs. How can I fix this issue? I also attached the screenshots. fft, frequency MATLAB Answers — New Questions
How so I use the Reinforcement Learning toolbox to train a Soft robot
Can I use the Reinforcement Learning toolbox to train a soft robot? None of the tutorials are for anything but rigid bots. Im not even sure if i can get a soft robot model into simulink as I couldnt figure it out earlierCan I use the Reinforcement Learning toolbox to train a soft robot? None of the tutorials are for anything but rigid bots. Im not even sure if i can get a soft robot model into simulink as I couldnt figure it out earlier Can I use the Reinforcement Learning toolbox to train a soft robot? None of the tutorials are for anything but rigid bots. Im not even sure if i can get a soft robot model into simulink as I couldnt figure it out earlier reinforcement learning, model, soft robot MATLAB Answers — New Questions
Can you use MATLAB with AutoCad?
I am interested in Architectural Acoustics, and I am wondering if you can use MATLAB with AutoCad in the design of recording studios and other spaces like concert venues, etc.?I am interested in Architectural Acoustics, and I am wondering if you can use MATLAB with AutoCad in the design of recording studios and other spaces like concert venues, etc.? I am interested in Architectural Acoustics, and I am wondering if you can use MATLAB with AutoCad in the design of recording studios and other spaces like concert venues, etc.? signal processing MATLAB Answers — New Questions
How to run code from a later line multiple times without restarting the whole script?
I have a script that is over a hundred lines long. For sake of argument, lets say lines 1 to 100 take 30 seconds to process.
I am currently working on lines 101 – 130. This 30 line block of code relies on the data that is processed from lines 1 to 100. I am new to MATLAB, so a lot of my coding from lines 101 – 130 encounters errors frequently.
The issue is I have to re-run the code from line 1 just to see if a new change I’ve made to lines 101 to 130 actually works.
I have tried using breakpoints by setting a breakpoint at line 100. But after making a change then continuing my script from line 100, I still have to run the code from line 1 if I encounter an error from line 100 to 130.
My question therefore, is can I somehow ‘cache’ the data from line 1 to 100? Is it therefore possible to run my code from line 101 to 130 multiple times over as I make constant changes through multiple errors without ever going back to line 1?I have a script that is over a hundred lines long. For sake of argument, lets say lines 1 to 100 take 30 seconds to process.
I am currently working on lines 101 – 130. This 30 line block of code relies on the data that is processed from lines 1 to 100. I am new to MATLAB, so a lot of my coding from lines 101 – 130 encounters errors frequently.
The issue is I have to re-run the code from line 1 just to see if a new change I’ve made to lines 101 to 130 actually works.
I have tried using breakpoints by setting a breakpoint at line 100. But after making a change then continuing my script from line 100, I still have to run the code from line 1 if I encounter an error from line 100 to 130.
My question therefore, is can I somehow ‘cache’ the data from line 1 to 100? Is it therefore possible to run my code from line 101 to 130 multiple times over as I make constant changes through multiple errors without ever going back to line 1? I have a script that is over a hundred lines long. For sake of argument, lets say lines 1 to 100 take 30 seconds to process.
I am currently working on lines 101 – 130. This 30 line block of code relies on the data that is processed from lines 1 to 100. I am new to MATLAB, so a lot of my coding from lines 101 – 130 encounters errors frequently.
The issue is I have to re-run the code from line 1 just to see if a new change I’ve made to lines 101 to 130 actually works.
I have tried using breakpoints by setting a breakpoint at line 100. But after making a change then continuing my script from line 100, I still have to run the code from line 1 if I encounter an error from line 100 to 130.
My question therefore, is can I somehow ‘cache’ the data from line 1 to 100? Is it therefore possible to run my code from line 101 to 130 multiple times over as I make constant changes through multiple errors without ever going back to line 1? matlab, breakpoints, cache MATLAB Answers — New Questions
