Error in MATLAB Code : BD18_Inverse Function

Question

Priyam Mishra 2023-10-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2032954-error-in-matlab-code-bd18_inverse-function

回答： SAI SRUJAN 2023-10-30

I am currently working on a MATLAB project involving the BD18_Inverse function for geotechnical data analysis, specifically the inverse filtering procedure. I am facing an issue with my MATLAB code and would appreciate your expertise in troubleshooting the problem.

The code executes the BD18_Inverse function, but I encounter an error that prevents the code from running successfully. The error message I receive is as follows:

```

Array indices must be positive integers or logical values.

Error in BD18_Inverse/CalculateW1 (line 79)

z50 = 1 + 2 * (C2 .* (z50_ref - 1)) .* (1 - 1 ./ (1 + ((q0t) ./ qinv(abs(z_prime)))).^0.5);

Error in BD18_Inverse (line 12)

w1 = CalculateW1(z_prime, n_m, dcone);

Error in BD18_Script (line 13)

BD18_Inverse(z, conetip, fsm, dcone);

```

Complete Code is here

function [qinv, fsinv] = BD18_Inverse(z, qm, fsm, dcone, sigma_v, Pa)
    % Initialize parameters
    n_m = 1;
    z_prime = (z - z(1)) / dcone;
    mt = 0.1;
    % Initialize qinv as qm
    qinv = qm;
    
    % Part 1: Inversion for Tip Resistance
    while true
        w1 = CalculateW1(z_prime, n_m, dcone);
        w2 = CalculateW2(qinv, n_m);
        wc = (w1 .* w2) / sum(w1 .* w2);
        
        % Calculate qtz' for smoothing
        qtz0 = mean(qinv(abs(z_prime) <= 0.05));
        for i = 1:length(z_prime)
            qtz0 = mean(qinv(abs(z_prime - z_prime(i)) <= 0.05));
            qz0t = sum(qinv .* wc);
            qinv(i) = qm(i) + (qinv(i) - qz0t) + (qtz0 - qz0t);
        end
        
        % Smoothing
        qinv = smooth(qinv, 5);
        
        % Calculate err
        err = sum(abs(qm - qinv)) / sum(qm);
        
        if err < 1e-6
            break;
        end
    end
    
    % Smoothing to remove high-frequency noise
    qinv = smooth(qinv, 5);
    
    % Part 2: Inversion for Sleeve Friction
    for i = 1:length(z)
        Q_tn_m = CalculateQt(qm, sigma_v, Pa, n_m);
        F_r_m = CalculateFr(fsm, qm, sigma_v);
        I_c_m = CalculateIc(Q_tn_m, F_r_m);
        n_m = CalculateN(I_c_m, sigma_v, Pa);
        % Check for convergence (n_m does not change significantly)
        if abs(n_m - 1) < 1e-6
            break;
        end
    end
    % Calculate Qtn for final iteration
    Q_tn_m = CalculateQt(qm, sigma_v, Pa, n_m);
    % Part 3: Interface Correction
    sharp_transitions = abs(log(qinv(2:end)) - log(qinv(1:end-1))) / (z_prime(2:end) - z_prime(1:end-1)) > mt;
    transition_zones = IdentifyTransitionZones(sharp_transitions, dcone, z_prime);
    qinv = PerformInterfaceCorrection(qinv, transition_zones);
    
    % Calculate fsinv
    Q_tn_inv = CalculateQt(qinv, sigma_v, Pa, n_m);
    F_r_inv = CalculateFr(fsm, qinv, sigma_v);
    I_c_inv = CalculateIc(Q_tn_inv, F_r_inv);
    n_inv = CalculateN(I_c_inv, sigma_v, Pa);
    
    fsinv = CalculateFsInv(F_r_inv, qinv, sigma_v);
    % Helper functions for Part 2
  function w1 = CalculateW1(z_prime, n_m, dcone)
    disp(size(z_prime)); % Add this line to check the size of z_prime
    disp(size(qinv));   % Add this line to check the size of qinv
    % Modify the code to match the dimensions
    z_prime = reshape(z_prime, size(qinv));
    z50_ref = 4.2;
    C2 = ones(size(z_prime));
    C2(z_prime >= 0) = 0;
    C2(z_prime >= 4) = 0.8;
    q0t = qinv(abs(z_prime) <= 0.05);
    z50 = 1 + (2 * (C2 .* (z50_ref - 1)) ).* (1 - 1 ./ (1 + ((q0t) ./ qinv(abs(z_prime)))).^0.5);
    w1 = 1 ./ (1 + (z_prime ./ z50).^3);
    w1(abs(z_prime) >= 4) = 0.5;
end
    function w2 = CalculateW2(qinv, n_m)
        qz0t = sum(qinv);
        w2 = sqrt(2 ./ (1 + (qz0t ./ qinv).^n_m));
    end
    function Q_tn = CalculateQt(q, sigma_v, Pa, n)
        Q_tn = ((q - sigma_v) ./ Pa) .* (Pa ./ sigma_v).^n;
    end
    function F_r = CalculateFr(fs, q, sigma_v)
        F_r = (fs ./ (q - sigma_v)) * 100;
    end
    function I_c = CalculateIc(Q_tn, F_r)
        I_c = sqrt((3.47 - log(Q_tn)).^2 + (log(F_r) + 1.22).^2);
    end
    function n = CalculateN(I_c, sigma_v, Pa)
        n = 0.381 * I_c + 0.05 * sigma_v / Pa - 0.15;
        n(n > 1) = 1;
    end
    % Helper functions for Part 3
    function transition_zones = IdentifyTransitionZones(sharp_transitions, dcone, z_prime)
        transition_zones = [];
        i = 1;
        while i <= length(sharp_transitions)
            if sharp_transitions(i) == 1
                transition_start = i;
                while i < length(sharp_transitions) && sharp_transitions(i) == 1
                    i = i + 1;
                end
                transition_end = i;
                transition_thickness = z_prime(transition_end) - z_prime(transition_start);
                if transition_thickness >= 3 * dcone
                    transition_zones = [transition_zones; transition_start, transition_end, transition_thickness];
                end
            end
            i = i + 1;
        end
    end
    function q = PerformInterfaceCorrection(q, transition_zones)
        for i = 1:size(transition_zones, 1)
            start_index = transition_zones(i, 1);
            end_index = transition_zones(i, 2);
            if q(end_index) > q(start_index)
                upper_zone = start_index + round(0.4 * transition_zones(i, 3));
                q(start_index:upper_zone) = q(upper_zone);
                
                lower_zone = upper_zone + 1;
                q(lower_zone:end_index) = q(end_index);
            else
                upper_zone = start_index + round(0.6 * transition_zones(i, 3));
                q(start_index:upper_zone) = q(upper_zone);
                lower_zone = upper_zone + 1;
                q(lower_zone:end_index) = q(end_index);
            end
        end
    end
end

I have tried several troubleshooting steps, including checking the variable sizes, confirming inputs, and debugging with smaller data subsets. However, I have not been able to identify the root cause of the issue.

To assist you in helping me, I have posted the MATLAB code in the following MATLAB Answers thread:

I kindly request your guidance and expertise in resolving this issue. Any insights, suggestions, or solutions you can provide would be greatly appreciated.

Thank you for your time and assistance.

Best regards,

Priyam Mishra

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

请先登录，再进行评论。

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

Answer 1

SAI SRUJAN 2023-10-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2032954-error-in-matlab-code-bd18_inverse-function#answer_1343436

在 MATLAB Online 中打开

Hi Priyam Mishra,

I understand that you are facing an issue involving the "BD18_Inverse" function for geotechnical data analysis, specifically the inverse filtering procedure.

The error message suggests that there is an issue with array indexing in the "CalculateW1" function.

The error is being caused by the line of code provided below.

z50 = 1 + 2 * (C2 .* (z50_ref - 1)) .* (1 - 1 ./ (1 + ((q0t) ./ qinv(abs(z_prime)))).^0.5);

The variable "z_prime" is utilized for indexing the "qinv" vector. It is essential to ensure that the value of "z_prime" adheres to the conditions, which require array indices to be positive integers or logical values. Additionally, it is advised to verify the absence of any potential division by zero error.

You can refer to the following documentation to understand more about "array indexing" in MATLAB,

https://in.mathworks.com/help/matlab/math/array-indexing.html

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

请先登录，再进行评论。

Error in MATLAB Code : BD18_Inverse Function

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

回答（1 个）

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

另请参阅

类别

标签

Community Treasure Hunt

Error in MATLAB Code : BD18_Inverse Function

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

回答（1 个）

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

另请参阅

类别

标签

Community Treasure Hunt

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

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