Help with Contour plot

i have this function in matlab

function [STR_SDV] = SIS_OP_SDV (K_SDV,T1_SDV)
     %variables are K_SDV and T1_SDV
     LS_SDV = 3.94E-6; 
     DC_SDV = 0.2;
     B_SDV = 0.02;
     MTTRSD_SDV =8;
     N_SDV=4;
     BS_SDV=B_SDV/2;
     DCS_SDV= LS_SDV*DC_SDV;
     LSU_SDV = (1-DCS_SDV)*LS_SDV ;
     LSD_SDV=DCS_SDV*LS_SDV ;
     C =(factorial (abs(N_SDV)))/(factorial(abs(K_SDV)-1))*((1-BS_SDV)*LSU_SDV + (1-BS_SDV)*LSD_SDV)^(K_SDV);
     l=1:(K_SDV-1);
     MDT_SDV_IND =(LSU_SDV/(LSU_SDV+LSD_SDV))*((T1_SDV./(l+1))+MTTRSD_SDV)+(LSD_SDV/(LSU_SDV+LSD_SDV))*MTTRSD_SDV; 
     PMDT_SDV_IND = prod(MDT_SDV_IND);
     STR_SDV_IND=C*PMDT_SDV_IND;
     STR_SDV_CCF=BS_SDV*LSU_SDV+BS_SDV*LSD_SDV;
     STR_SDV =STR_SDV_IND+STR_SDV_CCF;
 end

and i want to make a contour for it whereas the x are K_SDV and the y are T1_SDV And the K_SDV is between (1;3) and T1_SDV (4380;17520)

Answer 1

[Walter Roberson]

Some improvements in your code:

function [STR_SDV] = SIS_OP_SDV (K_SDV,T1_SDV)
     %variables are K_SDV and T1_SDV
     LS_SDV = 3.94E-6; 
     DC_SDV = 0.2;
     B_SDV = 0.02;
     MTTRSD_SDV = 8;
     N_SDV = 4;
     BS_SDV = B_SDV/2;
     DCS_SDV = LS_SDV.*DC_SDV;
     LSU_SDV = (1-DCS_SDV).*LS_SDV ;
     LSD_SDV = DCS_SDV.*LS_SDV ;
     C = (factorial (abs(N_SDV)))./(factorial(abs(K_SDV)-1)).*((1-BS_SDV).*LSU_SDV + (1-BS_SDV).*LSD_SDV).^(K_SDV);
     l = 1:(K_SDV-1);
     MDT_SDV_IND = (LSU_SDV./(LSU_SDV+LSD_SDV)).*((T1_SDV(:)./(l+1))+MTTRSD_SDV)+(LSD_SDV./(LSU_SDV+LSD_SDV)).*MTTRSD_SDV; 
     PMDT_SDV_IND = prod(MDT_SDV_IND,2);
     STR_SDV_IND = C.*PMDT_SDV_IND;
     STR_SDV_CCF = BS_SDV.*LSU_SDV+BS_SDV.*LSD_SDV;
     STR_SDV = STR_SDV_IND+STR_SDV_CCF;
     STD_SDV = reshape(STR_SDV, size(T1_SDV));
 end

With these changes, you must still call the function with a scalar K_SDV, but you can call with a non-scalar T1_SDV, and the result will be the same size as that T1_SDV.

Your code must use integer K_SDV -- you use factorial() and you use 1:K_SDV-1 both of which expect integer K_SDV. The 1:K_SDV-1 part impedes using non-scalar K_SDV, and since you use different array sizes for different K_SDV values it would be a bit of a nuisance to try to vectorize on K_SDV.

I emphasize again: you must use integer K_SDV. No linspace(1,3) for example.

T1_SDV does not have that restriction: you can use non-integer values for it if you want. linspace(4380, 17520) would be fine for that parameter.

Putting this all together, you can create a vector of T1_SDV values first, and then arrayfun over your integer K_SDV values.

It will be a poor contour indeed as you will only have 3 points along one of the axes.

You might consider using gamma(abs(K_SDV)-1+1) to allow non-integer K_SDV, but you still have the problem that l = 1:(K_SDV-1) which means you want a different number of terms in the prod() depending on K_SDV, and there just isn't any practical way to multiply a fractional number of terms if K_SDV were non-integer.

Reminder: 1:0 is empty so when K_SDV is 1, 1:(K_SDV-1) is 1:0 which is empty, leading to dividing T1_SDV by emptiness (since empty + 1 is empty), leading to MDT_SDV_IND being empty, and then prod(empty) is 1. So when K_SDV is 1, PMDT_SDV_IND will come out as 1. When K_SDV is 2 then 1:(K_SDV-1) would be 1:1 --> 1 so MDT_SDV_IND would have only one term for the product. The prod() does not start to become meaningful until K_SDV = 3 leading to 1:(3-1) --> 1:2

[Houcem eddine]

Thank you so much for your improvements I'm wondering if you can help me in making the contour plot for this function

[Walter Roberson]

As I indicated,

Putting this all together, you can create a vector of T1_SDV values first, and then arrayfun over your integer K_SDV values.

... put the results into a matrix, and contour() the matrix.

But as I indicated earlier, it will loop pretty poor unless you increase your maximum integer K_SDV value.

[Walter Roberson]

K_SDV_values = 1:3;
num_K_SDV = length(K_SDV_values);

T1_SDV_values = linspace(4380, 17520, 150);
num_T1_SDV = length(num_T1_SDV);

STR_SDVs = zeros(num_T1_SDV, num_K_SDV);  %pre-allocate

for K_SDV_idx = 1 : num_K_SDV
    K_SDV = K_SDV_values(K_SDV_idx);
    this_STR_SDV = SIS_OP_SDV(K_SDV,T1_SDV);
    STR_SDVs(:, K_SDV_idx) = this_STR_SDV;
end

surf(K_SDV_values, T1_SDV_values, STR_SDVs);