Category: News
NXP S32k1xx Toolbox installation
Hi everyone,
I’m trying to install nxp toolbox but getting given below error kindly help me out.
Dot indexing is not supported for variables of this type.
Error in NXP_Support_Package_S32K1xx>ConvertToolbox_Callback (line 143)
user_selection = get(handles.mbdt_version,’Value’);
^^^^^^^^^^^^^^^^^^^^
Error in gui_mainfcn (line 95)
feval(varargin{:});
^^^^^^^^^^^^^^^^^^
Error in NXP_Support_Package_S32K1xx (line 34)
gui_mainfcn(gui_State, varargin{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in matlab.graphics.internal.figfile.FigFile/read>@(hObject,eventdata)NXP_Support_Package_S32K1xx(‘ConvertToolbox_Callback’,hObject,eventdata,guidata(hObject))
Error using matlab.ui.internal.controller.uicontrol.UIControlController/triggerActionEvent (line 76)
Error while evaluating UIControl Callback.
regards,
Durgaprasad BHi everyone,
I’m trying to install nxp toolbox but getting given below error kindly help me out.
Dot indexing is not supported for variables of this type.
Error in NXP_Support_Package_S32K1xx>ConvertToolbox_Callback (line 143)
user_selection = get(handles.mbdt_version,’Value’);
^^^^^^^^^^^^^^^^^^^^
Error in gui_mainfcn (line 95)
feval(varargin{:});
^^^^^^^^^^^^^^^^^^
Error in NXP_Support_Package_S32K1xx (line 34)
gui_mainfcn(gui_State, varargin{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in matlab.graphics.internal.figfile.FigFile/read>@(hObject,eventdata)NXP_Support_Package_S32K1xx(‘ConvertToolbox_Callback’,hObject,eventdata,guidata(hObject))
Error using matlab.ui.internal.controller.uicontrol.UIControlController/triggerActionEvent (line 76)
Error while evaluating UIControl Callback.
regards,
Durgaprasad B Hi everyone,
I’m trying to install nxp toolbox but getting given below error kindly help me out.
Dot indexing is not supported for variables of this type.
Error in NXP_Support_Package_S32K1xx>ConvertToolbox_Callback (line 143)
user_selection = get(handles.mbdt_version,’Value’);
^^^^^^^^^^^^^^^^^^^^
Error in gui_mainfcn (line 95)
feval(varargin{:});
^^^^^^^^^^^^^^^^^^
Error in NXP_Support_Package_S32K1xx (line 34)
gui_mainfcn(gui_State, varargin{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in matlab.graphics.internal.figfile.FigFile/read>@(hObject,eventdata)NXP_Support_Package_S32K1xx(‘ConvertToolbox_Callback’,hObject,eventdata,guidata(hObject))
Error using matlab.ui.internal.controller.uicontrol.UIControlController/triggerActionEvent (line 76)
Error while evaluating UIControl Callback.
regards,
Durgaprasad B nxp toolbox MATLAB Answers — New Questions
How can I edit the properties for multiple requirements in the Requirements Editor in MATLAB R2024b?
Consider the requirement set file shown in the figure below, which is open in MATLAB R2024b. When I select multiple requirements, the "Properties" pane is greyed out and I cannot edit the properties for multiple requirements. How can I get around this?Consider the requirement set file shown in the figure below, which is open in MATLAB R2024b. When I select multiple requirements, the "Properties" pane is greyed out and I cannot edit the properties for multiple requirements. How can I get around this? Consider the requirement set file shown in the figure below, which is open in MATLAB R2024b. When I select multiple requirements, the "Properties" pane is greyed out and I cannot edit the properties for multiple requirements. How can I get around this? simulink, requirements, requirements-toolbox MATLAB Answers — New Questions
Help with Plotting the Envelope of a Scatter Plot in MATLAB
Hello everyone,
I’m trying to plot the envelope of a scatter plot in MATLAB. However, when I use the envelope function, it returns the same plot without highlighting the upper or lower bounds.
This plot was created using a scatter plot, and when I attempt to filter it to retain only the extremum values, I end up with just a single slice of the data.
Has anyone encountered a similar issue or knows how I could extract and plot the envelope correctly (in red), as shown in the reference image?
Thanks in advance!Hello everyone,
I’m trying to plot the envelope of a scatter plot in MATLAB. However, when I use the envelope function, it returns the same plot without highlighting the upper or lower bounds.
This plot was created using a scatter plot, and when I attempt to filter it to retain only the extremum values, I end up with just a single slice of the data.
Has anyone encountered a similar issue or knows how I could extract and plot the envelope correctly (in red), as shown in the reference image?
Thanks in advance! Hello everyone,
I’m trying to plot the envelope of a scatter plot in MATLAB. However, when I use the envelope function, it returns the same plot without highlighting the upper or lower bounds.
This plot was created using a scatter plot, and when I attempt to filter it to retain only the extremum values, I end up with just a single slice of the data.
Has anyone encountered a similar issue or knows how I could extract and plot the envelope correctly (in red), as shown in the reference image?
Thanks in advance! scatter, enveloppe, filtering MATLAB Answers — New Questions
I AM TRYING TO CONTROL SYNCHRONOUS RELUCTANCE MOTOR USING FOC. (Rs) = 8 Ω, (Ld) = 1.2 H (Lq) = 0.1 H P= 2, (J) = 0.125 kg·m² , (B) = 0.009 IF THE SIMULATION DIAGRAM IS OK
Post Content Post Content synchronous reluctance motor using foc MATLAB Answers — New Questions
How Microsoft Graph PowerShell SDK Access Tokens Work
Automatic Management of Access Tokens
Some years ago, I wrote about the access (bearer) tokens used by Entra ID. At the time, I focused on the access tokens obtained by apps from https://login.microsoftonline.com rather than those used by the Microsoft Graph PowerShell SDK.
One of the big advantages of using the Microsoft Graph PowerShell SDK is that developers don’t need to manage token renewal. When a script or app runs the Connect-MgGraph cmdlet to authenticate, an access token is obtained to allow cmdlets to run. When that access token approaches its expiration time, the Graph SDK requests a new token automatically.
Unless you knew that automatic renewal happens, you probably won’t realize how the Graph PowerShell SDK acquires and manages access tokens because details of the access token aren’t surfaced by a cmdlet like Get-MgContext. Although Get-MgContext reveals details of the current authentication context such as whether delegated or app-only authentication was used and the scopes (permissions) available to the session, there’s no trace of the access token.
Finding the Access Token Used by a Microsoft Graph PowerShell SDK Interactive Session
Some might be surprised that it’s not easier to find what access token is being used during a Graph PowerShell SDK session. However, automatic token management means that knowing what an access token is and when the token will expire is not information that’s necessary for a session to function, so it’s reasonable to keep the data hidden behind the scenes.
To find the access token, it’s necessary to make a special form of request to any Graph API. The request can be made in an interactive session or an app-only session. This example uses a request against the drives endpoint to retrieve retention label information for a file, but any request will work:
$Uri = ("https://graph.microsoft.com/v1.0/drives/{0}/items/{1}/retentionLabel" -f $OneDriveInfo.Id, $File.Id)
$Data = Invoke-MgGraphRequest -Uri $Uri -Method Get -OutputType HttpResponseMessage
The key point here is that the Invoke-MgGraphRequest cmdlet specifies that it should receive a HTTP response rather than data. The request specified in the URI is simply a way to ask the Graph for the HTTP response. The response contains several interesting components:
Version : 2.0
Content : System.Net.Http.DecompressionHandler+GZipDecompressedContent
StatusCode : OK
ReasonPhrase : OK
Headers : {[Cache-Control, System.String[]], [Vary, System.String[]], [Strict-Transport-Security, System.String[]], [request-id, System.String[]]…}
TrailingHeaders : {}
RequestMessage : Method: GET, RequestUri: 'https://graph.microsoft.com/v1.0/drives/b!_xwZzApnQEeEWOYGdTfHR_FlEFWmBHl JixksigwWMZ_hpEW05Pd_R7OzPT4YdqXq/items/01R343MZ43HNLCSCCT3ZBLLUIJGB3GJ5B3/retentionLabel',
Version: 2.0, Content: <null>, Headers:
{
User-Agent: Mozilla/5.0
User-Agent: (Windows NT 10.0; Microsoft Windows 10.0.26100; en-IE)
User-Agent: PowerShell/7.5.2
User-Agent: Invoke-MgGraphRequest
FeatureFlag: 00000003
Cache-Control: no-store, no-cache
Authorization: Bearer eyJ0eXAiOiJKV1QiLCJub25jZSI6IlZrTmh0QjdFajZpSUhRVkRwdmZYeVVldUEyeFFBbFhyR1M
The access token is at the bottom of the output and can be retrieved with:
$Data.RequestMessage.Headers.Authorization.Parameter
Decrypting an Access Token
Isolating the access token makes it easier to copy and input into the jwt.io token decoder. Figure 1 shows the raw JSON output; selecting the claims table tab presents the information in a more easily understood fashion (this reference page also helps).
The decoded token reveals details like the app in use (Microsoft Graph Command Line Tools), its identifier, the user, and the available permissions
"app_displayname": "Microsoft Graph Command Line Tools", "appid": "14d82eec-204b-4c2f-b7e8-296a70dab67e", "family_name": "Redmond", "given_name": "Tony", "idtyp": "user", "ipaddr": "109.78.233.203", "name": "Tony Redmond", "oid": "eff4cd58-1bb8-4899-94de-795f656b4a18", "scp": "AccessReview.Read.All Agreement.Read.All Analytics.Read APIConnectors.Read.All Application.Read.All Application.ReadWrite.All AppRoleAssignment.ReadWrite.All AuditLog.Read.All AuditLogsQuery.Read.All ...
The token also contain timestamps in UNIX epoch format for when the token was issued and when it will expire. The claims table output shows the date in local time. You can also convert these dates with PowerShell:
$UnixEpochValue = 1752763429 $Date = [DateTimeOffset]::FromUnixTimeSeconds($UnixEpochValue).ToLocalTime().DateTime Write-Host "UNIX epoch $UnixEpochValue is" $(Get-Date $Date -format 'dd-MMM-yyyy HH:mm') UNIX epoch 1752763429 is 17-Jul-2025 15:43
Down further in the token you’ll find the wids array, which holds the identifiers for the Entra ID roles held by the user. Remember, during an interactive Graph SDK session the available permissions are the intersection between delegated permissions and administrative roles. In other words, if access to data isn’t available through a permission, it might be through a role.
Reusing a Graph Access Token
You can take the access token used by the Graph interactive session and use it to retrieve information without using a Graph SDK cmdlet. In this code snippet, we prepare a hash table containing the access token formatted in the way that Graph requests expect the data to be presented and use the token with the Invoke-RestMethod cmdlet to find the details of the signed-in user.
$Headers = @{}
$Headers.Add("Authorization", ("{0} {1}" -f $Data.RequestMessage.Headers.Authorization.Scheme, $Data.RequestMessage.Headers.Authorization.Parameter))
$Me = Invoke-RestMethod -Uri "https://graph.microsoft.com/v1.0/me" -Headers $Headers
$Me.displayName
Tony Redmond
Interesting But Not Very Useful Information
All of this is firmly in the interesting but not very useful category. If an app wants to make Graph API requests without using the Microsoft Graph PowerShell SDK, it will do the norm and obtain an access token programmatically before running any requests. The permissions available to the app are the set of delegated and application permissions held by the app’s service principal. If the app runs for over an hour, it will need to renew the access token.
Apart from testing code to write this article, I don’t think I have ever looked at the access token in a Microsoft Graph PowerShell SDK session. I might in the future, but right now I can’t think of a good reason why I should.
Need some assistance to write and manage PowerShell scripts for Microsoft 365? 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.
Error when applying coder to WLAN Toolbox example code
An error occurs when using Matlab coder with hWLANPacketDetector.m used in SDRBeaconReceiverExample.mlx.An error occurs when using Matlab coder with hWLANPacketDetector.m used in SDRBeaconReceiverExample.mlx. An error occurs when using Matlab coder with hWLANPacketDetector.m used in SDRBeaconReceiverExample.mlx. wlan coder MATLAB Answers — New Questions
Please guide me, I can’t find the F280049 Flash file, on TMS320F280049 when using Uniflash it will only work once, after reset it won’t work.
Please guide me, I can’t find the F280049 Flash file, on TMS320F280049 when using Uniflash it will only work once, after reset it won’t work. After reset, it might be from the MCU boosting from RAM?Please guide me, I can’t find the F280049 Flash file, on TMS320F280049 when using Uniflash it will only work once, after reset it won’t work. After reset, it might be from the MCU boosting from RAM? Please guide me, I can’t find the F280049 Flash file, on TMS320F280049 when using Uniflash it will only work once, after reset it won’t work. After reset, it might be from the MCU boosting from RAM? c2000, simulink MATLAB Answers — New Questions
Expansion of expressions involving Einstein summation convention with Levi-Civita tensors
Dear all
I am evaluating expressions of the type
where all the subindices , and the dot on the variable , , symbolizes the time derivative.
I’d be interested in being able to ask MatLab for the fully expanded , , and , so I can then numerically evaluate those terms. Ideally, it would be a plus if it also gave me the monstrous explicit expressions in LaTeX and/or could extrapolate the expanded expression to C++.
Do you have any ideas on how I could do this?Dear all
I am evaluating expressions of the type
where all the subindices , and the dot on the variable , , symbolizes the time derivative.
I’d be interested in being able to ask MatLab for the fully expanded , , and , so I can then numerically evaluate those terms. Ideally, it would be a plus if it also gave me the monstrous explicit expressions in LaTeX and/or could extrapolate the expanded expression to C++.
Do you have any ideas on how I could do this? Dear all
I am evaluating expressions of the type
where all the subindices , and the dot on the variable , , symbolizes the time derivative.
I’d be interested in being able to ask MatLab for the fully expanded , , and , so I can then numerically evaluate those terms. Ideally, it would be a plus if it also gave me the monstrous explicit expressions in LaTeX and/or could extrapolate the expanded expression to C++.
Do you have any ideas on how I could do this? symbolic to numeric expressions, einstein summation convention MATLAB Answers — New Questions
Resonance in moonpool type structures
I have the following code for my problem where i apply the match eigen function technique to find the unknown coefficients and then plotting one result wave motion inside a moonpool structure. However the expected plot for me should be decreasing after the resonance peak. But in my code it is decrease but then again blows up. I did not find any source of error. The analytical expressions which i used for my problem is write in pdf file at the end. Here is my code . I just want to figure out my source of error to fix it but i couldn’t. So I need help in acheiving my desired results. I also attached the refrence plot.
clear all;
close all;
clc;
tic;
% Parameters from your setup
N = 5; % Number of evanescent modes
M = 3; % Number of azimuthal modes
RE = 1.0; % Outer radius (m)
ratio = 2.0; % Ratio R_E / R_I
RI = RE / ratio; % Inner radius (m)
h = 4.0 * RE; % Water depth (m)
b = 2.0; % Distance from cylinder bottom to seabed (m)
d = h – b; % Interface depth (m)
g = 9.81; % Gravity (m/s^2)
l = 0; % Azimuthal order
% I_given_wavenumber = 1;
X_kRE = linspace(0.01, 4, 200);
% Initialize eta0 array before loop
eta0 = zeros(size(X_kRE)); % Preallocate
for idx = 1:length(X_kRE)
% if I_given_wavenumber == 1
wk = X_kRE(idx) / RE; % dimensional wavenumber
omega = sqrt(g * wk * tanh(wk * h)); % dispersion relation
fun_alpha = @(x) omega^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes
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);
end
end
% % find evanescent modes k1…k_{N-1} stored in k_all(2)…k_all(N)
% for n = 2 : N
% m = n – 1; % mode index for evanescent
% a = ((2*m-1)*pi)/(2*h) + eps;
% b = (m*pi)/h – eps;
% % at a: tan→ -∞ so f(a)<0; at b: tan→0 so f(b)=+omega^2>0
% k_all(n) = fzero(@(k) g*k.*tan(k*h) + omega^2, [a, b]);
% 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
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Define right-hand side vector G
G = zeros(2*M+2*N, 1); % Total length = M + N+ M+ N = 2M+2N
%
% Block 1: G (size M = 4)
for i = 1:M
if i == 1
G(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
G(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2))*besselj(l, wk*RE); %Z_m^l
end
end
%
% Block 2: Y (size N = 3)
for j = 1:N
if j == 1
G(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
G(j + M, 1) = 0; % Y_n^l
end
end
% Block 3: X (size M = 4)
for i = 1:M
G(i + M + N, 1) = 0; % X_0^l, X_m^l
end
% Block 4: W (size N = 3)
for j = 1:N
G(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);
% Construct the full block matrix H
H = [ I_M, -A, ZMM, ZMN; % Block row 1
-B, I_N, -C, ZNN; % Block row 2
ZMM, ZMN, I_M, -D; % Block row 3
-E, ZNN, -F, I_N]; % Block row 4
% X = H G; % Solve the linear system
X = pinv(H) * G; % Use pseudoinverse for stability
%
b_vec = X(1:M); % Coefficients b^l
a_vec = X(M+1 : M+N); % Coefficients a^l
c_vec = X(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = X(2*M+N+1 : end); % Coefficients d^l
%%%%%%%%% 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
eta0(idx) = abs(term1+sum1);
% % disp(k.’) % Should all be real and positive for n ≥ 2
% %
% %
end
plot(X_kRE, eta0, ‘k’, ‘LineWidth’, 1.5);
xlabel(‘$k_0 R_E$’, ‘Interpreter’, ‘latex’);
ylabel(‘$eta / (iA)$’, ‘Interpreter’, ‘latex’);
title(‘Wave motion amplitude for RE/RI = 4.0’);
grid on;
xlim([0 4]); % Match the reference plot
ylim([0 6]); % Optional: based on expected peak height
elapsedTime = toc;
disp([‘Time consuming = ‘, num2str(elapsedTime), ‘ s’]);
function out = dbesselh(l, z)
if l == 0
out = -besselh(1, 1, z);
else
out = 0.5 * (besselh(l – 1, 1, z) – besselh(l + 1, 1, z));
end
end
function out = dbesseli(l, z)
if l == 0
out = besseli(1, z);
else
out = 0.5 * (besseli(l – 1, z) + besseli(l + 1, z));
end
end
function out = dbesselj(l, z)
if l == 0
out = -besselj(1, z);
else
out = 0.5 * (besselj(l – 1, z) – besselj(l + 1, z));
end
end
function out = dbesselk(l, z)
if l == 0
out = -besselk(1, z);
else
out = -0.5 * (besselk(l – 1, z) + besselk(l + 1, z));
end
endI have the following code for my problem where i apply the match eigen function technique to find the unknown coefficients and then plotting one result wave motion inside a moonpool structure. However the expected plot for me should be decreasing after the resonance peak. But in my code it is decrease but then again blows up. I did not find any source of error. The analytical expressions which i used for my problem is write in pdf file at the end. Here is my code . I just want to figure out my source of error to fix it but i couldn’t. So I need help in acheiving my desired results. I also attached the refrence plot.
clear all;
close all;
clc;
tic;
% Parameters from your setup
N = 5; % Number of evanescent modes
M = 3; % Number of azimuthal modes
RE = 1.0; % Outer radius (m)
ratio = 2.0; % Ratio R_E / R_I
RI = RE / ratio; % Inner radius (m)
h = 4.0 * RE; % Water depth (m)
b = 2.0; % Distance from cylinder bottom to seabed (m)
d = h – b; % Interface depth (m)
g = 9.81; % Gravity (m/s^2)
l = 0; % Azimuthal order
% I_given_wavenumber = 1;
X_kRE = linspace(0.01, 4, 200);
% Initialize eta0 array before loop
eta0 = zeros(size(X_kRE)); % Preallocate
for idx = 1:length(X_kRE)
% if I_given_wavenumber == 1
wk = X_kRE(idx) / RE; % dimensional wavenumber
omega = sqrt(g * wk * tanh(wk * h)); % dispersion relation
fun_alpha = @(x) omega^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes
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);
end
end
% % find evanescent modes k1…k_{N-1} stored in k_all(2)…k_all(N)
% for n = 2 : N
% m = n – 1; % mode index for evanescent
% a = ((2*m-1)*pi)/(2*h) + eps;
% b = (m*pi)/h – eps;
% % at a: tan→ -∞ so f(a)<0; at b: tan→0 so f(b)=+omega^2>0
% k_all(n) = fzero(@(k) g*k.*tan(k*h) + omega^2, [a, b]);
% 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
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Define right-hand side vector G
G = zeros(2*M+2*N, 1); % Total length = M + N+ M+ N = 2M+2N
%
% Block 1: G (size M = 4)
for i = 1:M
if i == 1
G(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
G(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2))*besselj(l, wk*RE); %Z_m^l
end
end
%
% Block 2: Y (size N = 3)
for j = 1:N
if j == 1
G(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
G(j + M, 1) = 0; % Y_n^l
end
end
% Block 3: X (size M = 4)
for i = 1:M
G(i + M + N, 1) = 0; % X_0^l, X_m^l
end
% Block 4: W (size N = 3)
for j = 1:N
G(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);
% Construct the full block matrix H
H = [ I_M, -A, ZMM, ZMN; % Block row 1
-B, I_N, -C, ZNN; % Block row 2
ZMM, ZMN, I_M, -D; % Block row 3
-E, ZNN, -F, I_N]; % Block row 4
% X = H G; % Solve the linear system
X = pinv(H) * G; % Use pseudoinverse for stability
%
b_vec = X(1:M); % Coefficients b^l
a_vec = X(M+1 : M+N); % Coefficients a^l
c_vec = X(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = X(2*M+N+1 : end); % Coefficients d^l
%%%%%%%%% 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
eta0(idx) = abs(term1+sum1);
% % disp(k.’) % Should all be real and positive for n ≥ 2
% %
% %
end
plot(X_kRE, eta0, ‘k’, ‘LineWidth’, 1.5);
xlabel(‘$k_0 R_E$’, ‘Interpreter’, ‘latex’);
ylabel(‘$eta / (iA)$’, ‘Interpreter’, ‘latex’);
title(‘Wave motion amplitude for RE/RI = 4.0’);
grid on;
xlim([0 4]); % Match the reference plot
ylim([0 6]); % Optional: based on expected peak height
elapsedTime = toc;
disp([‘Time consuming = ‘, num2str(elapsedTime), ‘ s’]);
function out = dbesselh(l, z)
if l == 0
out = -besselh(1, 1, z);
else
out = 0.5 * (besselh(l – 1, 1, z) – besselh(l + 1, 1, z));
end
end
function out = dbesseli(l, z)
if l == 0
out = besseli(1, z);
else
out = 0.5 * (besseli(l – 1, z) + besseli(l + 1, z));
end
end
function out = dbesselj(l, z)
if l == 0
out = -besselj(1, z);
else
out = 0.5 * (besselj(l – 1, z) – besselj(l + 1, z));
end
end
function out = dbesselk(l, z)
if l == 0
out = -besselk(1, z);
else
out = -0.5 * (besselk(l – 1, z) + besselk(l + 1, z));
end
end I have the following code for my problem where i apply the match eigen function technique to find the unknown coefficients and then plotting one result wave motion inside a moonpool structure. However the expected plot for me should be decreasing after the resonance peak. But in my code it is decrease but then again blows up. I did not find any source of error. The analytical expressions which i used for my problem is write in pdf file at the end. Here is my code . I just want to figure out my source of error to fix it but i couldn’t. So I need help in acheiving my desired results. I also attached the refrence plot.
clear all;
close all;
clc;
tic;
% Parameters from your setup
N = 5; % Number of evanescent modes
M = 3; % Number of azimuthal modes
RE = 1.0; % Outer radius (m)
ratio = 2.0; % Ratio R_E / R_I
RI = RE / ratio; % Inner radius (m)
h = 4.0 * RE; % Water depth (m)
b = 2.0; % Distance from cylinder bottom to seabed (m)
d = h – b; % Interface depth (m)
g = 9.81; % Gravity (m/s^2)
l = 0; % Azimuthal order
% I_given_wavenumber = 1;
X_kRE = linspace(0.01, 4, 200);
% Initialize eta0 array before loop
eta0 = zeros(size(X_kRE)); % Preallocate
for idx = 1:length(X_kRE)
% if I_given_wavenumber == 1
wk = X_kRE(idx) / RE; % dimensional wavenumber
omega = sqrt(g * wk * tanh(wk * h)); % dispersion relation
fun_alpha = @(x) omega^2 + g*x.*tan(x*h); %dispersion equation for evanescent modes
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);
end
end
% % find evanescent modes k1…k_{N-1} stored in k_all(2)…k_all(N)
% for n = 2 : N
% m = n – 1; % mode index for evanescent
% a = ((2*m-1)*pi)/(2*h) + eps;
% b = (m*pi)/h – eps;
% % at a: tan→ -∞ so f(a)<0; at b: tan→0 so f(b)=+omega^2>0
% k_all(n) = fzero(@(k) g*k.*tan(k*h) + omega^2, [a, b]);
% 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
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Define right-hand side vector G
G = zeros(2*M+2*N, 1); % Total length = M + N+ M+ N = 2M+2N
%
% Block 1: G (size M = 4)
for i = 1:M
if i == 1
G(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ; % Z_0^l
else
G(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2))*besselj(l, wk*RE); %Z_m^l
end
end
%
% Block 2: Y (size N = 3)
for j = 1:N
if j == 1
G(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2); % Y_0^l
else
G(j + M, 1) = 0; % Y_n^l
end
end
% Block 3: X (size M = 4)
for i = 1:M
G(i + M + N, 1) = 0; % X_0^l, X_m^l
end
% Block 4: W (size N = 3)
for j = 1:N
G(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);
% Construct the full block matrix H
H = [ I_M, -A, ZMM, ZMN; % Block row 1
-B, I_N, -C, ZNN; % Block row 2
ZMM, ZMN, I_M, -D; % Block row 3
-E, ZNN, -F, I_N]; % Block row 4
% X = H G; % Solve the linear system
X = pinv(H) * G; % Use pseudoinverse for stability
%
b_vec = X(1:M); % Coefficients b^l
a_vec = X(M+1 : M+N); % Coefficients a^l
c_vec = X(M+N+1 : 2*M+N); % Coefficients c^l
d_vec = X(2*M+N+1 : end); % Coefficients d^l
%%%%%%%%% 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
eta0(idx) = abs(term1+sum1);
% % disp(k.’) % Should all be real and positive for n ≥ 2
% %
% %
end
plot(X_kRE, eta0, ‘k’, ‘LineWidth’, 1.5);
xlabel(‘$k_0 R_E$’, ‘Interpreter’, ‘latex’);
ylabel(‘$eta / (iA)$’, ‘Interpreter’, ‘latex’);
title(‘Wave motion amplitude for RE/RI = 4.0’);
grid on;
xlim([0 4]); % Match the reference plot
ylim([0 6]); % Optional: based on expected peak height
elapsedTime = toc;
disp([‘Time consuming = ‘, num2str(elapsedTime), ‘ s’]);
function out = dbesselh(l, z)
if l == 0
out = -besselh(1, 1, z);
else
out = 0.5 * (besselh(l – 1, 1, z) – besselh(l + 1, 1, z));
end
end
function out = dbesseli(l, z)
if l == 0
out = besseli(1, z);
else
out = 0.5 * (besseli(l – 1, z) + besseli(l + 1, z));
end
end
function out = dbesselj(l, z)
if l == 0
out = -besselj(1, z);
else
out = 0.5 * (besselj(l – 1, z) – besselj(l + 1, z));
end
end
function out = dbesselk(l, z)
if l == 0
out = -besselk(1, z);
else
out = -0.5 * (besselk(l – 1, z) + besselk(l + 1, z));
end
end moonpool, resonance, system of linear equation, semi analytical solution, wave motion MATLAB Answers — New Questions
How we can calculculate the MTF from PSF?
I want to calculate the MTF from this PSF.I want to calculate the MTF from this PSF. I want to calculate the MTF from this PSF. frequency MATLAB Answers — New Questions
The provided solution is coming up as incorrect.
I am currently doing the MATLAB onramp course and I am on task 6 of the final project Stellar moon 2 and I cant figure out the answer. When I select the solution I copied it and put it as my answer however it comes up as incorrect. This is beyond frusterating as I am unable to progress and I have tried deleting everything and trying again several times and I just want to finish this course. Can someone please help.I am currently doing the MATLAB onramp course and I am on task 6 of the final project Stellar moon 2 and I cant figure out the answer. When I select the solution I copied it and put it as my answer however it comes up as incorrect. This is beyond frusterating as I am unable to progress and I have tried deleting everything and trying again several times and I just want to finish this course. Can someone please help. I am currently doing the MATLAB onramp course and I am on task 6 of the final project Stellar moon 2 and I cant figure out the answer. When I select the solution I copied it and put it as my answer however it comes up as incorrect. This is beyond frusterating as I am unable to progress and I have tried deleting everything and trying again several times and I just want to finish this course. Can someone please help. solve, onramp MATLAB Answers — New Questions
speeding up the data analysis
I am working on a project where I have to analyse high frequency wind data (50 Hz) coming from a wind turbine data measurement system. this amount of data must be converted from .dat binary files to .mat files which I can use in matlab. The data must then be filtered and then averaged over 10 minutes to be compared to the data from another measurement system. Doing all this requires the analysis of thousands of data and it’s very time consuming (right now it is about 105 seconds just for the data of 2 days). How can I speed the process up?
for ii = 3:length(filedir)
filename = filedir(ii).name;
newdata = ReadFamosDataIntoTimeTable(filename);
% filter
IsConsidered = newdata.avbladeangleGRe<40 … % normal operation
& newdata.RAWS9>0.5 … % good availability
& ~isnan(newdata.iv10mswindspeed2GRe) & newdata.av100msabswinddirectionGRe>180 & newdata.av100msabswinddirectionGRe<250;
TurbineData(ii-2).date = filename;
TurbineData(ii-2).iv10mswindspeed2GRe = mean(newdata.iv10mswindspeed2GRe(IsConsidered));
TurbineData(ii-2).ivactivepowerGRe = mean(newdata.ivactivepowerGRe(IsConsidered));
TurbineData(ii-2).CalculatedAirdensity_GRe = mean(newdata.CalculatedAirdensity_GRe(IsConsidered));
TurbineData(ii-2).HWShub1 = mean(newdata.HWShub1(IsConsidered));
TurbineData(ii-2).av100msabswinddirectionGRe = mean(newdata.av100msabswinddirectionGRe(IsConsidered));
end
The function ReadFamosDataIntoTimeTable is to convert the .dat binary data into .mat dataI am working on a project where I have to analyse high frequency wind data (50 Hz) coming from a wind turbine data measurement system. this amount of data must be converted from .dat binary files to .mat files which I can use in matlab. The data must then be filtered and then averaged over 10 minutes to be compared to the data from another measurement system. Doing all this requires the analysis of thousands of data and it’s very time consuming (right now it is about 105 seconds just for the data of 2 days). How can I speed the process up?
for ii = 3:length(filedir)
filename = filedir(ii).name;
newdata = ReadFamosDataIntoTimeTable(filename);
% filter
IsConsidered = newdata.avbladeangleGRe<40 … % normal operation
& newdata.RAWS9>0.5 … % good availability
& ~isnan(newdata.iv10mswindspeed2GRe) & newdata.av100msabswinddirectionGRe>180 & newdata.av100msabswinddirectionGRe<250;
TurbineData(ii-2).date = filename;
TurbineData(ii-2).iv10mswindspeed2GRe = mean(newdata.iv10mswindspeed2GRe(IsConsidered));
TurbineData(ii-2).ivactivepowerGRe = mean(newdata.ivactivepowerGRe(IsConsidered));
TurbineData(ii-2).CalculatedAirdensity_GRe = mean(newdata.CalculatedAirdensity_GRe(IsConsidered));
TurbineData(ii-2).HWShub1 = mean(newdata.HWShub1(IsConsidered));
TurbineData(ii-2).av100msabswinddirectionGRe = mean(newdata.av100msabswinddirectionGRe(IsConsidered));
end
The function ReadFamosDataIntoTimeTable is to convert the .dat binary data into .mat data I am working on a project where I have to analyse high frequency wind data (50 Hz) coming from a wind turbine data measurement system. this amount of data must be converted from .dat binary files to .mat files which I can use in matlab. The data must then be filtered and then averaged over 10 minutes to be compared to the data from another measurement system. Doing all this requires the analysis of thousands of data and it’s very time consuming (right now it is about 105 seconds just for the data of 2 days). How can I speed the process up?
for ii = 3:length(filedir)
filename = filedir(ii).name;
newdata = ReadFamosDataIntoTimeTable(filename);
% filter
IsConsidered = newdata.avbladeangleGRe<40 … % normal operation
& newdata.RAWS9>0.5 … % good availability
& ~isnan(newdata.iv10mswindspeed2GRe) & newdata.av100msabswinddirectionGRe>180 & newdata.av100msabswinddirectionGRe<250;
TurbineData(ii-2).date = filename;
TurbineData(ii-2).iv10mswindspeed2GRe = mean(newdata.iv10mswindspeed2GRe(IsConsidered));
TurbineData(ii-2).ivactivepowerGRe = mean(newdata.ivactivepowerGRe(IsConsidered));
TurbineData(ii-2).CalculatedAirdensity_GRe = mean(newdata.CalculatedAirdensity_GRe(IsConsidered));
TurbineData(ii-2).HWShub1 = mean(newdata.HWShub1(IsConsidered));
TurbineData(ii-2).av100msabswinddirectionGRe = mean(newdata.av100msabswinddirectionGRe(IsConsidered));
end
The function ReadFamosDataIntoTimeTable is to convert the .dat binary data into .mat data optimization, for loop, speed, data conversion MATLAB Answers — New Questions
how to buy student lisence i am in Egypt and all credit card i tried not works.
i want to buy student lisence, is there any restrictions on
Egypt country? all credit card i tried not works at your web site.i want to buy student lisence, is there any restrictions on
Egypt country? all credit card i tried not works at your web site. i want to buy student lisence, is there any restrictions on
Egypt country? all credit card i tried not works at your web site. specify libraries, or specific apis, matlab, license MATLAB Answers — New Questions
Switching values around in a matrix
Hi,
Say I have a 5 by 2 matrix in the form:
A = [2 5; 9 7; 10 2; 3 2; 1 9]
And I want to make a 10 by 1 matrix from these values so that I get:
B = [2;5;9;7;10;2;3;2;1;9]
How would I do this?
I know there is a probably a simple fix, but I haven’t been able to do it.
Many thanks,
ScottHi,
Say I have a 5 by 2 matrix in the form:
A = [2 5; 9 7; 10 2; 3 2; 1 9]
And I want to make a 10 by 1 matrix from these values so that I get:
B = [2;5;9;7;10;2;3;2;1;9]
How would I do this?
I know there is a probably a simple fix, but I haven’t been able to do it.
Many thanks,
Scott Hi,
Say I have a 5 by 2 matrix in the form:
A = [2 5; 9 7; 10 2; 3 2; 1 9]
And I want to make a 10 by 1 matrix from these values so that I get:
B = [2;5;9;7;10;2;3;2;1;9]
How would I do this?
I know there is a probably a simple fix, but I haven’t been able to do it.
Many thanks,
Scott matrices, matrix manipulation MATLAB Answers — New Questions
Why do I receive the error “Unable to save login information. You are currently signed in but you will be prompted to sign in again the next time that you start this application”?
When launching MATLAB, why do I receive the error "Unable to save login information. You are currently signed in but you will be prompted to sign in again the next time that you start this application"?When launching MATLAB, why do I receive the error "Unable to save login information. You are currently signed in but you will be prompted to sign in again the next time that you start this application"? When launching MATLAB, why do I receive the error "Unable to save login information. You are currently signed in but you will be prompted to sign in again the next time that you start this application"? MATLAB Answers — New Questions
Plot showing gaps in numerical solution
I’m working on a numerical model involving the merging of four turbulent plumes, each originating from a circular area source located at the corners of a square. My goal is to compute and visualize the merged outer contour that encloses all four interacting plumes. I have attached the code below, the equation I am solving is at the bottom of this code. When I make the plot some gaps appear and is there a better way to plot it so the gaps disappear
%% Four-Plume Contour + Area Plotting with Automatic k_crit Detection
close all; clc; clf
global r0 k theta
%%Parameter Input
r0 = 0.25; %Plume Source Radius
%% STEP 1: Find k_crit by setting θ = π/4 and solving corresponding function, g(rho) = 0 for r0
theta = pi/4;
kL = 0*(1 – r0^2)^2; % Conservative lower bound
kU = 1.05*(1 – r0^2)^2; % Previous Theoretical maximum
k_testArray = linspace(kL,kU, 501);
rho_valid = NaN(size(k_testArray)); % Store valid roots for each k in sweep
k_crit = NaN; % First k where g(rho)=0 has a real root
contactDetected = false; % Logical flag to stop at first contact
for jj = 1:length(k_testArray)
k = k_testArray(jj);
rho = mnhf_secant(@FourPlumes_Equation, [0.1 0.5], 1e-6, 0);
if rho > 0 && isreal(rho) && isfinite(rho)
rho_valid(jj) = rho;
if ~contactDetected
k_crit = k;
contactDetected = true;
end
end
end
fprintf(‘Detected critical k value (first contact with θ = π/4): k_crit = %.8fnn’, k_crit);
%% Parameter input: k values
Nt = 10001; % Resolution, must be an odd number
theta_array = linspace(-pi/2, pi/2, Nt);
k_array = [0.6]
%k_array = [0.2 0.3 0.4 0.5 0.6 0.7 k_crit 0.8 0.9 1.2 1.5]; % test range values for k
area = zeros(1,length(k_array));
%% STEP 2: Loop over k values, plot, and compute area
figure(1); hold on; box on
for jj = 1:length(k_array)
k = k_array(jj);
% For merged Case
if k >= k_crit
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square; xlim([0 2]); ylim([-1 1]);
% Plot θ = π/4 lines for visual reference
x_diag = linspace(0, 2, Nt); y_diag = x_diag;
plot(x_diag, y_diag, ‘r–‘); plot(x_diag, -y_diag, ‘r–‘);
% Area estimation using trapezoidal rule
iq = find(y1 >= y_diag); % intersection
if ~isempty(iq)
area(jj) = 2 * abs(trapz(x1(min(iq):(Nt-1)/2+1), y1(min(iq):(Nt-1)/2+1)));
area(jj) = area(jj) + x1(min(iq)-1)^2;
end
else
% Before Merging Case
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square;
% Detect breakdown region and recompute in "gap"
ij = find(isnan(rho_array1(1:round(end/2))));
if ~isempty(ij)
theta_min = theta_array(max(ij));
theta_array2 = linspace(theta_min, -theta_min, Nt);
rho_array2 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array2(ii);
rho_array2(ii) = mnhf_secant(@FourPlumes_Equation, [0.1 0.9], 1e-6, 0);
if rho_array2(ii) < 0
rho_array2(ii) = NaN;
end
end
[x2, y2] = pol2cart(theta_array2, rho_array2);
plot(x2, y2, ‘r.’);
ik = (Nt-1)/2 + find(isnan(x2((Nt-1)/2+1:end)));
x = [x2((Nt-1)/2+1:min(ik)-1), x1(max(ij)+1:(Nt-1)/2+1)];
y = [y2((Nt-1)/2+1:min(ik)-1), -y1(max(ij)+1:(Nt-1)/2+1)];
plot(x, y, ‘b–‘); plot(x, -y, ‘b–‘);
% Area from symmetric region
area(jj) = 2 * trapz(x, y);
end
end
% Plot plume source circles (n = 4)
pos1 = [1 – r0, -r0, 2*r0, 2*r0]; % (1, 0)
pos2 = [-1 – r0, -r0, 2*r0, 2*r0]; % (-1, 0)
pos3 = [-r0, 1 – r0, 2*r0, 2*r0]; % (0, 1)
pos4 = [-r0, -1 – r0, 2*r0, 2*r0]; % (0, -1)
rectangle(‘Position’, pos1, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos2, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos3, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos4, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
end
%% Function: Four-Plume Equation (Li & Flynn B3) —
function f = FourPlumes_Equation(rho)
global r0 k theta
f = (rho.^2 – 2*rho*cos(theta) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi/2) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + 1.5*pi) + 1 – r0^2) – k^2;
end
function r=mnhf_secant(Fun,x,tol,trace)
%MNHF_SECANT finds the root of "Fun" using secant scheme.
%
% Fun – name of the external function
% x – vector of length 2, (initial guesses)
% tol – tolerance
% trace – print intermediate results
%
% Usage mnhf_secant(@poly1,[-0.5 2.0],1e-8,1)
% poly1 is the name of the external function.
% [-0.5 2.0] are the initial guesses for the root.
% Check inputs.
if nargin < 4, trace = 1; end
if nargin < 3, tol = 1e-8; end
if (length(x) ~= 2)
error(‘Please provide two initial guesses’)
end
f = feval(Fun,x); % Fun is assumed to accept a vector
for i = 1:100
x3 = x(1)-f(1)*(x(2)-x(1))/(f(2)-f(1)); % Update the guess.
f3 = feval(Fun,x3); % Function evaluation.
if trace, fprintf(1,’%3i %12.5f %12.5fn’, i,x3,f3); end
if abs(f3) < tol % Check for convergenece.
r = x3;
return
else % Reset values for x(1), x(2), f(1) and f(2).
x(1) = x(2); f(1) = f(2); x(2) = x3; f(2) = f3;
end
end
r = NaN;
endI’m working on a numerical model involving the merging of four turbulent plumes, each originating from a circular area source located at the corners of a square. My goal is to compute and visualize the merged outer contour that encloses all four interacting plumes. I have attached the code below, the equation I am solving is at the bottom of this code. When I make the plot some gaps appear and is there a better way to plot it so the gaps disappear
%% Four-Plume Contour + Area Plotting with Automatic k_crit Detection
close all; clc; clf
global r0 k theta
%%Parameter Input
r0 = 0.25; %Plume Source Radius
%% STEP 1: Find k_crit by setting θ = π/4 and solving corresponding function, g(rho) = 0 for r0
theta = pi/4;
kL = 0*(1 – r0^2)^2; % Conservative lower bound
kU = 1.05*(1 – r0^2)^2; % Previous Theoretical maximum
k_testArray = linspace(kL,kU, 501);
rho_valid = NaN(size(k_testArray)); % Store valid roots for each k in sweep
k_crit = NaN; % First k where g(rho)=0 has a real root
contactDetected = false; % Logical flag to stop at first contact
for jj = 1:length(k_testArray)
k = k_testArray(jj);
rho = mnhf_secant(@FourPlumes_Equation, [0.1 0.5], 1e-6, 0);
if rho > 0 && isreal(rho) && isfinite(rho)
rho_valid(jj) = rho;
if ~contactDetected
k_crit = k;
contactDetected = true;
end
end
end
fprintf(‘Detected critical k value (first contact with θ = π/4): k_crit = %.8fnn’, k_crit);
%% Parameter input: k values
Nt = 10001; % Resolution, must be an odd number
theta_array = linspace(-pi/2, pi/2, Nt);
k_array = [0.6]
%k_array = [0.2 0.3 0.4 0.5 0.6 0.7 k_crit 0.8 0.9 1.2 1.5]; % test range values for k
area = zeros(1,length(k_array));
%% STEP 2: Loop over k values, plot, and compute area
figure(1); hold on; box on
for jj = 1:length(k_array)
k = k_array(jj);
% For merged Case
if k >= k_crit
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square; xlim([0 2]); ylim([-1 1]);
% Plot θ = π/4 lines for visual reference
x_diag = linspace(0, 2, Nt); y_diag = x_diag;
plot(x_diag, y_diag, ‘r–‘); plot(x_diag, -y_diag, ‘r–‘);
% Area estimation using trapezoidal rule
iq = find(y1 >= y_diag); % intersection
if ~isempty(iq)
area(jj) = 2 * abs(trapz(x1(min(iq):(Nt-1)/2+1), y1(min(iq):(Nt-1)/2+1)));
area(jj) = area(jj) + x1(min(iq)-1)^2;
end
else
% Before Merging Case
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square;
% Detect breakdown region and recompute in "gap"
ij = find(isnan(rho_array1(1:round(end/2))));
if ~isempty(ij)
theta_min = theta_array(max(ij));
theta_array2 = linspace(theta_min, -theta_min, Nt);
rho_array2 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array2(ii);
rho_array2(ii) = mnhf_secant(@FourPlumes_Equation, [0.1 0.9], 1e-6, 0);
if rho_array2(ii) < 0
rho_array2(ii) = NaN;
end
end
[x2, y2] = pol2cart(theta_array2, rho_array2);
plot(x2, y2, ‘r.’);
ik = (Nt-1)/2 + find(isnan(x2((Nt-1)/2+1:end)));
x = [x2((Nt-1)/2+1:min(ik)-1), x1(max(ij)+1:(Nt-1)/2+1)];
y = [y2((Nt-1)/2+1:min(ik)-1), -y1(max(ij)+1:(Nt-1)/2+1)];
plot(x, y, ‘b–‘); plot(x, -y, ‘b–‘);
% Area from symmetric region
area(jj) = 2 * trapz(x, y);
end
end
% Plot plume source circles (n = 4)
pos1 = [1 – r0, -r0, 2*r0, 2*r0]; % (1, 0)
pos2 = [-1 – r0, -r0, 2*r0, 2*r0]; % (-1, 0)
pos3 = [-r0, 1 – r0, 2*r0, 2*r0]; % (0, 1)
pos4 = [-r0, -1 – r0, 2*r0, 2*r0]; % (0, -1)
rectangle(‘Position’, pos1, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos2, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos3, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos4, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
end
%% Function: Four-Plume Equation (Li & Flynn B3) —
function f = FourPlumes_Equation(rho)
global r0 k theta
f = (rho.^2 – 2*rho*cos(theta) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi/2) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + 1.5*pi) + 1 – r0^2) – k^2;
end
function r=mnhf_secant(Fun,x,tol,trace)
%MNHF_SECANT finds the root of "Fun" using secant scheme.
%
% Fun – name of the external function
% x – vector of length 2, (initial guesses)
% tol – tolerance
% trace – print intermediate results
%
% Usage mnhf_secant(@poly1,[-0.5 2.0],1e-8,1)
% poly1 is the name of the external function.
% [-0.5 2.0] are the initial guesses for the root.
% Check inputs.
if nargin < 4, trace = 1; end
if nargin < 3, tol = 1e-8; end
if (length(x) ~= 2)
error(‘Please provide two initial guesses’)
end
f = feval(Fun,x); % Fun is assumed to accept a vector
for i = 1:100
x3 = x(1)-f(1)*(x(2)-x(1))/(f(2)-f(1)); % Update the guess.
f3 = feval(Fun,x3); % Function evaluation.
if trace, fprintf(1,’%3i %12.5f %12.5fn’, i,x3,f3); end
if abs(f3) < tol % Check for convergenece.
r = x3;
return
else % Reset values for x(1), x(2), f(1) and f(2).
x(1) = x(2); f(1) = f(2); x(2) = x3; f(2) = f3;
end
end
r = NaN;
end I’m working on a numerical model involving the merging of four turbulent plumes, each originating from a circular area source located at the corners of a square. My goal is to compute and visualize the merged outer contour that encloses all four interacting plumes. I have attached the code below, the equation I am solving is at the bottom of this code. When I make the plot some gaps appear and is there a better way to plot it so the gaps disappear
%% Four-Plume Contour + Area Plotting with Automatic k_crit Detection
close all; clc; clf
global r0 k theta
%%Parameter Input
r0 = 0.25; %Plume Source Radius
%% STEP 1: Find k_crit by setting θ = π/4 and solving corresponding function, g(rho) = 0 for r0
theta = pi/4;
kL = 0*(1 – r0^2)^2; % Conservative lower bound
kU = 1.05*(1 – r0^2)^2; % Previous Theoretical maximum
k_testArray = linspace(kL,kU, 501);
rho_valid = NaN(size(k_testArray)); % Store valid roots for each k in sweep
k_crit = NaN; % First k where g(rho)=0 has a real root
contactDetected = false; % Logical flag to stop at first contact
for jj = 1:length(k_testArray)
k = k_testArray(jj);
rho = mnhf_secant(@FourPlumes_Equation, [0.1 0.5], 1e-6, 0);
if rho > 0 && isreal(rho) && isfinite(rho)
rho_valid(jj) = rho;
if ~contactDetected
k_crit = k;
contactDetected = true;
end
end
end
fprintf(‘Detected critical k value (first contact with θ = π/4): k_crit = %.8fnn’, k_crit);
%% Parameter input: k values
Nt = 10001; % Resolution, must be an odd number
theta_array = linspace(-pi/2, pi/2, Nt);
k_array = [0.6]
%k_array = [0.2 0.3 0.4 0.5 0.6 0.7 k_crit 0.8 0.9 1.2 1.5]; % test range values for k
area = zeros(1,length(k_array));
%% STEP 2: Loop over k values, plot, and compute area
figure(1); hold on; box on
for jj = 1:length(k_array)
k = k_array(jj);
% For merged Case
if k >= k_crit
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square; xlim([0 2]); ylim([-1 1]);
% Plot θ = π/4 lines for visual reference
x_diag = linspace(0, 2, Nt); y_diag = x_diag;
plot(x_diag, y_diag, ‘r–‘); plot(x_diag, -y_diag, ‘r–‘);
% Area estimation using trapezoidal rule
iq = find(y1 >= y_diag); % intersection
if ~isempty(iq)
area(jj) = 2 * abs(trapz(x1(min(iq):(Nt-1)/2+1), y1(min(iq):(Nt-1)/2+1)));
area(jj) = area(jj) + x1(min(iq)-1)^2;
end
else
% Before Merging Case
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
if rho_array1(ii) < 0
rho_array1(ii) = NaN;
end
end
[x1, y1] = pol2cart(theta_array, rho_array1);
plot(x1, y1, ‘-k’); axis square;
% Detect breakdown region and recompute in "gap"
ij = find(isnan(rho_array1(1:round(end/2))));
if ~isempty(ij)
theta_min = theta_array(max(ij));
theta_array2 = linspace(theta_min, -theta_min, Nt);
rho_array2 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array2(ii);
rho_array2(ii) = mnhf_secant(@FourPlumes_Equation, [0.1 0.9], 1e-6, 0);
if rho_array2(ii) < 0
rho_array2(ii) = NaN;
end
end
[x2, y2] = pol2cart(theta_array2, rho_array2);
plot(x2, y2, ‘r.’);
ik = (Nt-1)/2 + find(isnan(x2((Nt-1)/2+1:end)));
x = [x2((Nt-1)/2+1:min(ik)-1), x1(max(ij)+1:(Nt-1)/2+1)];
y = [y2((Nt-1)/2+1:min(ik)-1), -y1(max(ij)+1:(Nt-1)/2+1)];
plot(x, y, ‘b–‘); plot(x, -y, ‘b–‘);
% Area from symmetric region
area(jj) = 2 * trapz(x, y);
end
end
% Plot plume source circles (n = 4)
pos1 = [1 – r0, -r0, 2*r0, 2*r0]; % (1, 0)
pos2 = [-1 – r0, -r0, 2*r0, 2*r0]; % (-1, 0)
pos3 = [-r0, 1 – r0, 2*r0, 2*r0]; % (0, 1)
pos4 = [-r0, -1 – r0, 2*r0, 2*r0]; % (0, -1)
rectangle(‘Position’, pos1, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos2, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos3, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos4, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
end
%% Function: Four-Plume Equation (Li & Flynn B3) —
function f = FourPlumes_Equation(rho)
global r0 k theta
f = (rho.^2 – 2*rho*cos(theta) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi/2) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + 1.5*pi) + 1 – r0^2) – k^2;
end
function r=mnhf_secant(Fun,x,tol,trace)
%MNHF_SECANT finds the root of "Fun" using secant scheme.
%
% Fun – name of the external function
% x – vector of length 2, (initial guesses)
% tol – tolerance
% trace – print intermediate results
%
% Usage mnhf_secant(@poly1,[-0.5 2.0],1e-8,1)
% poly1 is the name of the external function.
% [-0.5 2.0] are the initial guesses for the root.
% Check inputs.
if nargin < 4, trace = 1; end
if nargin < 3, tol = 1e-8; end
if (length(x) ~= 2)
error(‘Please provide two initial guesses’)
end
f = feval(Fun,x); % Fun is assumed to accept a vector
for i = 1:100
x3 = x(1)-f(1)*(x(2)-x(1))/(f(2)-f(1)); % Update the guess.
f3 = feval(Fun,x3); % Function evaluation.
if trace, fprintf(1,’%3i %12.5f %12.5fn’, i,x3,f3); end
if abs(f3) < tol % Check for convergenece.
r = x3;
return
else % Reset values for x(1), x(2), f(1) and f(2).
x(1) = x(2); f(1) = f(2); x(2) = x3; f(2) = f3;
end
end
r = NaN;
end numerical solution, nonlinear-equations, root-finding, polar-coordinates, contour-gaps MATLAB Answers — New Questions
Monthly Update #122 Available for Office 365 for IT Pros eBook
August 2025 Update Available for Subscribers to Download
The Office 365 for IT Pros team is delighted to announce the availability of monthly update #122 for Office 365 for IT Pros (2026 edition). This is the first update since the publication of the 2026 book, and it includes a bunch of changes covering new information unearthed or researched over the last month. You can see more details about the updates in our change log.
Subscribers for the 2026 edition can fetch the updated files using the link in the receipt emailed to them from Gumroad.com after their purchase. The link is personalized to the subscriber and always fetches the latest files. See our FAQ for more information.
Many previous subscribers have renewed for the 2026 edition. We are humbled by the level of support we receive from the technical community. If you are a previous subscriber, you should have received email from us with instructions about how to extend your subscription at a heavily discounted rate. Because we became tired of Gumroad.com emails being blocked as spam, we sent individual emails to thousands of previous subscribers. I believe every subscriber should have received a note from us by now. If not, please contact us at o365itprosrenewals AT office365itpros.com.
Microsoft Cloud Keeps on Growing
Microsoft released their FY25 Q4 results on July 30. One of the eye-popping figures was a big leap in Microsoft Cloud revenue to $46.7 billion (up 27% year-over-year). That’s $21 billion more in a quarter than Microsoft earned in FY23 Q1 and represents an annualized run rate of $186.8 billion. Microsoft 365 continues to grow strongly with its revenue up 16% year-over-year. Seat growth is slowing and was 6% year-over-year, mostly in SME and frontline workers. Slower seat growth is inevitable given the size of the installed base (over 430 million, according to CFO Amy Hood when discussing the Microsoft FY25 Q3 results). The growth in revenues is due to upsell to more expensive licenses, including Microsoft 365 Copilot.
Speaking about Copilot, in the best traditional of misleading figures given out during results briefings, Satya Nadella said “Our family of Copilot apps has surpassed 100 million monthly active users across commercial and consumer.” The number claimed sounds impressive, but what everyone really wants to know is how many Microsoft 365 Copilot paid seats are in active use. It was followed by the assertion that “Purview is used by three quarters of Microsoft 365 Copilot customers to protect their data.” Again, no real data to measure anything against.
Interestingly, once again Microsoft didn’t give an updated number for Teams users. The number remains at the 320 million monthly active users claimed in October 2023.
Speaking of numbers, Microsoft reported a 99.995% performance against the Microsoft 365 SLA for the second quarter of 2025. That might come as news to any tenant that experienced a significant outage between April and June, but the result is unsurprising given the number of tenants and seats. It takes a massive outage involving tens of millions of seats over several hours to budge the SLA needle slightly. Outages happen all the time, but none at the level of severity necessary to impact the SLA
Update for the Automating Microsoft 365 with PowerShell eBook
As is our practice, we released an update for the Automating Microsoft 365 with PowerShell eBook earlier than the monthly Office 365 for IT Pros update. The PowerShell book is included with Office 365 for IT Pros, so subscribers should see 4 files when they access the update on Gumroad.com (PDF and EPUB files for both books).
The PowerShell book is available separately, but Office 365 for IT Pros subscribers do not have to purchase it. If you want to read a physical copy, you can buy a paperback version from Amazon.com. I haven’t actually seen the paperback yet, but I plan to get some author copies for delivery to the TEC 2025 conference in Minneapolis at the end of September. Come to TEC and you might get a signed copy! If not, you can still attend TEC sessions delivered by Tony, Paul, Brian, and Michel.
On to Monthly Update #123
Microsoft 365 doesn’t stop changing and we don’t stop analyzing and documenting the most important technical information for Microsoft 365 tenant administrators. It’s what we’ve done since 2015.
Backslash (mldivide) slower than inverse and multiplication
The common wisdom is that Ay is more accurate than inv(A)*y and I believe that to be true, but when is it also faster? Say the matrix A is well-conditioned so I don’t really care about the ability of to find a least-squares solution.
Am I doing something misleading here? The takes much longer.
A = randn(20);
A = A*A’;
s = randn(20,1);
timeit(@() inv(A)*s)
timeit(@() A s)
From the documentation for mldivide, it sounds like it should be using the Cholesky solver since A is Hermitian, why is that not faster than an inv and matrix multiplication?
ishermitian(A)The common wisdom is that Ay is more accurate than inv(A)*y and I believe that to be true, but when is it also faster? Say the matrix A is well-conditioned so I don’t really care about the ability of to find a least-squares solution.
Am I doing something misleading here? The takes much longer.
A = randn(20);
A = A*A’;
s = randn(20,1);
timeit(@() inv(A)*s)
timeit(@() A s)
From the documentation for mldivide, it sounds like it should be using the Cholesky solver since A is Hermitian, why is that not faster than an inv and matrix multiplication?
ishermitian(A) The common wisdom is that Ay is more accurate than inv(A)*y and I believe that to be true, but when is it also faster? Say the matrix A is well-conditioned so I don’t really care about the ability of to find a least-squares solution.
Am I doing something misleading here? The takes much longer.
A = randn(20);
A = A*A’;
s = randn(20,1);
timeit(@() inv(A)*s)
timeit(@() A s)
From the documentation for mldivide, it sounds like it should be using the Cholesky solver since A is Hermitian, why is that not faster than an inv and matrix multiplication?
ishermitian(A) mldivide, inv, slow MATLAB Answers — New Questions
How to integrate PEM electrolysis system and PEM fuel cell system to form the DC Electrical Power System?
Fig. 1 shows the PEM (Polymer Electrolyte Membrane) Electrolysis System.
Fig.1 PEM Electrolysis System
Fig. 2 shows the PEM (Polymer Electrolyte Membrane) Fuel Cell System.
Fig.2 PEM Fuel Cell System
Fig. 3 shows the Solid-Oxide Fuel Cell Connected to Three-Phase Electrical Power System.
Fig.3 Solid-Oxide Fuel Cell Connected to Three-Phase Electrical Power System
My question is how to complete the following tasks:
(1) Integrate Fig. 1 and Fig. 2 to form the DC Electrical Power System—the PEM Electrolysis—Fuel Cell.
(2) Integrate Fig. 1, Fig. 2, and Fig. 3 to form the DC Electrical Power System—PEM Electrolysis—Fuel Cell—Electrical Power System.
Reference：
Hydrogen Electrolyzer: https://ww2.mathworks.cn/discovery/hydrogen-electrolyzer.html
PEM Fuel Cell System：https://ww2.mathworks.cn/help/simscape/ug/pem-fuel-cell-system.html
Fuel Cell Model：https://ww2.mathworks.cn/discovery/fuel-cell-model.html
Solid-Oxide Fuel Cell Connected to Three-Phase Electrical Power System：
https://ww2.mathworks.cn/help/sps/ug/solid-oxide-fuel-cell-connected-to-three-phase-electrical-power-system.html
Hydrogen Electrolyzer: https://ww2.mathworks.cn/discovery/hydrogen-electrolyzer.html
https://ww2.mathworks.cn/help/simscape/ug/pem-electrolysis-system.html
Integrating the PEM electrolysis system and the PEM fuel cell system—does this mean replacing the Hydrogen Source in the PEM fuel cell system with the PEM Electrolysis System? Could the experts provide specific guidance on this? Thanks in advance.Fig. 1 shows the PEM (Polymer Electrolyte Membrane) Electrolysis System.
Fig.1 PEM Electrolysis System
Fig. 2 shows the PEM (Polymer Electrolyte Membrane) Fuel Cell System.
Fig.2 PEM Fuel Cell System
Fig. 3 shows the Solid-Oxide Fuel Cell Connected to Three-Phase Electrical Power System.
Fig.3 Solid-Oxide Fuel Cell Connected to Three-Phase Electrical Power System
My question is how to complete the following tasks:
(1) Integrate Fig. 1 and Fig. 2 to form the DC Electrical Power System—the PEM Electrolysis—Fuel Cell.
(2) Integrate Fig. 1, Fig. 2, and Fig. 3 to form the DC Electrical Power System—PEM Electrolysis—Fuel Cell—Electrical Power System.
Reference：
Hydrogen Electrolyzer: https://ww2.mathworks.cn/discovery/hydrogen-electrolyzer.html
PEM Fuel Cell System：https://ww2.mathworks.cn/help/simscape/ug/pem-fuel-cell-system.html
Fuel Cell Model：https://ww2.mathworks.cn/discovery/fuel-cell-model.html
Solid-Oxide Fuel Cell Connected to Three-Phase Electrical Power System：
https://ww2.mathworks.cn/help/sps/ug/solid-oxide-fuel-cell-connected-to-three-phase-electrical-power-system.html
Hydrogen Electrolyzer: https://ww2.mathworks.cn/discovery/hydrogen-electrolyzer.html
https://ww2.mathworks.cn/help/simscape/ug/pem-electrolysis-system.html
Integrating the PEM electrolysis system and the PEM fuel cell system—does this mean replacing the Hydrogen Source in the PEM fuel cell system with the PEM Electrolysis System? Could the experts provide specific guidance on this? Thanks in advance. Fig. 1 shows the PEM (Polymer Electrolyte Membrane) Electrolysis System.
Fig.1 PEM Electrolysis System
Fig. 2 shows the PEM (Polymer Electrolyte Membrane) Fuel Cell System.
Fig.2 PEM Fuel Cell System
Fig. 3 shows the Solid-Oxide Fuel Cell Connected to Three-Phase Electrical Power System.
Fig.3 Solid-Oxide Fuel Cell Connected to Three-Phase Electrical Power System
My question is how to complete the following tasks:
(1) Integrate Fig. 1 and Fig. 2 to form the DC Electrical Power System—the PEM Electrolysis—Fuel Cell.
(2) Integrate Fig. 1, Fig. 2, and Fig. 3 to form the DC Electrical Power System—PEM Electrolysis—Fuel Cell—Electrical Power System.
Reference：
Hydrogen Electrolyzer: https://ww2.mathworks.cn/discovery/hydrogen-electrolyzer.html
PEM Fuel Cell System：https://ww2.mathworks.cn/help/simscape/ug/pem-fuel-cell-system.html
Fuel Cell Model：https://ww2.mathworks.cn/discovery/fuel-cell-model.html
Solid-Oxide Fuel Cell Connected to Three-Phase Electrical Power System：
https://ww2.mathworks.cn/help/sps/ug/solid-oxide-fuel-cell-connected-to-three-phase-electrical-power-system.html
Hydrogen Electrolyzer: https://ww2.mathworks.cn/discovery/hydrogen-electrolyzer.html
https://ww2.mathworks.cn/help/simscape/ug/pem-electrolysis-system.html
Integrating the PEM electrolysis system and the PEM fuel cell system—does this mean replacing the Hydrogen Source in the PEM fuel cell system with the PEM Electrolysis System? Could the experts provide specific guidance on this? Thanks in advance. pem electrolysis system, pem fuel cell system, solid-oxide fuel cell MATLAB Answers — New Questions
How to change the background color of uiradiobuttons?
I have a few uibuttongroups in my app and to distinguish them visually I want to apply different background colors to them. But the radiobuttons’ backgroundcolor can’t be changed so the whole setup looks a bit weird.
Here is a sample code:
fig=uifigure;
bg=uibuttongroup(fig,’BackgroundColor’,[0.8 0.9 1],’Position’,[100 100 150 150]);
rb1=uiradiobutton(bg,’Text’,’Previous’,’Position’,[10 100 75 20]);
rb2=uiradiobutton(bg,’Text’,’Next’,’Position’,[10 50 75 20]);
%rb1.BackgroundColor=[0.8 0.9 1]; %Won’t work!
In addition, the default background grey of uibuttongroup & uiradiobutton are slightly different in R2025a so even with the default colours, it feels like something is off!I have a few uibuttongroups in my app and to distinguish them visually I want to apply different background colors to them. But the radiobuttons’ backgroundcolor can’t be changed so the whole setup looks a bit weird.
Here is a sample code:
fig=uifigure;
bg=uibuttongroup(fig,’BackgroundColor’,[0.8 0.9 1],’Position’,[100 100 150 150]);
rb1=uiradiobutton(bg,’Text’,’Previous’,’Position’,[10 100 75 20]);
rb2=uiradiobutton(bg,’Text’,’Next’,’Position’,[10 50 75 20]);
%rb1.BackgroundColor=[0.8 0.9 1]; %Won’t work!
In addition, the default background grey of uibuttongroup & uiradiobutton are slightly different in R2025a so even with the default colours, it feels like something is off! I have a few uibuttongroups in my app and to distinguish them visually I want to apply different background colors to them. But the radiobuttons’ backgroundcolor can’t be changed so the whole setup looks a bit weird.
Here is a sample code:
fig=uifigure;
bg=uibuttongroup(fig,’BackgroundColor’,[0.8 0.9 1],’Position’,[100 100 150 150]);
rb1=uiradiobutton(bg,’Text’,’Previous’,’Position’,[10 100 75 20]);
rb2=uiradiobutton(bg,’Text’,’Next’,’Position’,[10 50 75 20]);
%rb1.BackgroundColor=[0.8 0.9 1]; %Won’t work!
In addition, the default background grey of uibuttongroup & uiradiobutton are slightly different in R2025a so even with the default colours, it feels like something is off! appdesigner, uiradiobuttons MATLAB Answers — New Questions
