Code an equation in MATLAB

In matlab, I write the following code to the above equation as follows

ri = sqrt(psi_coordinates.^2+theta_coordinates.^2);
A = sum((ri*delta_epsilon)/2);
V = delta_z*sum(A); 

But the result I obtained is not correct. Can anyone comment if I made some mistake please ?

Thank you .

[Walter Roberson]

you forgot the square root

[HN]

Thank you Walter, It was actually typo. I didn’t forget the square root in the code here. But the problem still existed. So, you think if the square root is included no other issues? Thank you again

[Dyuman Joshi]

I am assuming m is the number of rows and n is the number of columns.

If so, then you have to use sum((ri*delta_epsilon)/2,2) while defining A to get the row-wise sum.

[Jeffrey Clark]

Where is the sqrt in computation of ri?

ri = sqrt(psi_coordinates.^2+theta_coordinates.^2);

[Jeffrey Clark]

The last equation in the image is missing a closing paren in the first definition of V before simplification; I assume the open paren should be moved before the second summation and another close paren added before the simplification of V. V is supposed to be calculating a scaler sum of sums, which must be implying that the ri and A are to be working on only one of the dimensions of phi and theta, in MATLAB the sums of the ri and A will be only on one of them, and the V sum should then be on the other. The calculation of the sum in A in MATLAB will be on one of the dimensions. Since you didn't tell MATLAB which dimention the sum is on could you be calculating with the delta_epsilon and delta_z affecting the wrong dimentional data?

[HN]

Hello everyone. I am sorry for the delay and I really apperciate the supports. I was trying to follow your advices but nothing is changed. Then, I took some time to resolve it and no progress. Below I attached the data for psi and theta in text file. delta_epsilon is 0.0622. delta_z is also given as 0.0125. The expected result for V is 0.083794. I apperciate if someone check it and help me to figure out the mistake I made.

Just to make things more clear, the equations are

@Walter Roberson, @Dyuman Joshi @Jeffrey Clark

Thank you all

[Dyuman Joshi]

While defining delta_epsilon, what is the value of m?

[HN]

Hi @Dyuman Joshi,

value of m is 101 and number of n is 80. That is why the size of theta and psi is 80 x 101.

Thank you

[Dyuman Joshi]

Try this

r = (th.^2+psi.^2).^0.5;
delz=0.0125;
deleps=0.0622;

A = sum(r*deleps*delz/2,2); %row sum
V = delz*sum(A)

[HN]

@Dyuman Joshi, Thank you so much. I am running it but why is the squre root ignored?

Thanks again

[HN]

@Dyuman Joshi, the result I got with the above code, without including sqrt, is 7.9115e-04. The expected result is 0.083794. I am now trying it with sqrt while computing r. I really don't understand where I made a mistake.

Thanks

[Dyuman Joshi]

Sorry, it was ignored because of my mistake. I forgot to take it into account. I have edited it.

Did you take element wise sqrt? Like I have, in the edit?

[HN]

Yes, I take it into account but the result is far from the expected result. I also changed the element selection to column sum. But the result is not correct as well.

[Dyuman Joshi]

How are you inputting the data? I can try it as well and check.

[HN]

I saved the data attached in here from the workspace. When you load it, theta and psi will have the size of 80 by 101. Then use them to calculate V. That is how I used it. Thanks

[Dyuman Joshi]

@HN, I am getting 81*102 matrices while loading from both the text files.

[HN]

hmmm, How that happened? But the value of V is still far from the true value that was supposed to be.

[Dyuman Joshi]

Idk, but that's what I got. Even changing delesp according to m, brought only a minimum change in the value of V.

How do you get the value of V btw?

[HN]

I didn’t get it. It is from a paper.

Below is the part of program @Dyuman Joshi

while point(3) <= z_end + 0.000001
    work_points = zeros(2,1);
    z_vector = [z_vector point(3)];
    delta_meridian = 2*pi/meridian_number;
    for meridian = 0 : delta_meridian : 2*pi
        point(4) = 0;				
        point(5) = 0;				
        
        	if point(3) < 0.2*leg_length					
        		 delta_angular = delta_start / 200;			
        	elseif point(3) >= 0.2*leg_length && point(3) <= 0.3*leg_length	
        		 delta_angular =  delta_start / 50;		
        	else								
        		 delta_angular = delta_start;
            end
        
       %delta_angular = delta_start;		
        
        ccontinue = 1;				
        point(4) = point(4) + delta_angular * cos(meridian);
        point(5) = point(5) + delta_angular * sin(meridian);
        while  ccontinue == 1
            ccontinue = workspace(constrain(point),max_cond);
            if ccontinue == 0 && delta_angular > tolerance		
                point(4) = point(4) - delta_angular * cos(meridian);
                point(5) = point(5) - delta_angular * sin(meridian);
                delta_angular = delta_angular / 2;			
                ccontinue = 1;					
                
            end
            psi_temp = point(4); theta_temp = point(5);
            point(4) = point(4) + delta_angular * cos(meridian);
            point(5) = point(5) + delta_angular * sin(meridian);
        end
        temp_vect = [psi_temp;theta_temp];
        work_points = [work_points temp_vect];
    end
    work_surf = [work_surf' work_points']';
    point(3) = point(3) + delta_z;
    
    
    
end
size_temp = size(work_surf);
work_surf = work_surf(:,2:size_temp(2));
size_temp = size(work_surf);

psi_coordinates = []; theta_coordinates = [];
for i = 1 :2:(z_slices*2+2)			
    psi_coordinates = [psi_coordinates work_surf(i,:)'];
    theta_coordinates = [theta_coordinates work_surf(i+1,:)'];
end
psi_coordinates = psi_coordinates';
theta_coordinates = theta_coordinates';

 psi = psi_coordinates'; 
 th = theta_coordinates';
z_position = ones(z_slices+1,size_temp(2));	
for i=1:z_slices+1
    z_position(i,:) = z_vector(i);
end

ri = sqrt(psi.^2 + th.^2); 
A = sum((ri*delta_meridian*delta_z)/2,2);
v = delta_z*sum(A)

[Jeffrey Clark]

In your last code segment, next-to-last line references Delta_z not delta_z what is going on here? Retyping code incorrectly for these comments or somewhere have a capital D delta?

[HN]

Thank you @Jeffrey Clark,

It was typo when I write it here for those comments. There is no Delta_z and it is corrected.

Thank you

[Jeffrey Clark]

So the line below may not be doing what you expect due to rounding to computer precision; although there seems to be more missing than provided in your edited part of the program, this for loop will dictate the size of work_points and work_surf but I don't know to what extent the psi and th since z_slices value isn't given and may not jive with the generated work_surf:

for meridian = 0 : delta_meridian : 2*pi

[HN]

@Jeffrey Clark, Thank you. Your opinion gave me some clue, I guess.

Normally,ψ and θ have components,i.e., and . So, parameter r is obtained only from the terms. How do you think I can make a selection while calculating r ?

Many thanks in advance.

[Jeffrey Clark]

1. I think you could be helped most by uploading the complete MATLAB program
2. in your 10 May 2022 at 7:53 post, you added a delta_z in the calculation of A which wasn't in your original post
3. in your 10 May 2022 at 20:22 post, in the first if..elseif..else statement the else will never execute because the elseif includes the negated if clause and an or (||) with another clause that won't be tested - perhaps you intended an and (&&) in the elseif?
4. based on what you did provide, phi and th would not be matrices just vectors, so the sums don't make sense

The program seems too complex for something that that that just needs to create some sort of angle or sine space mesh.

[Dyuman Joshi]

What is the unit of meridian? Degrees/Radians?

size_temp = size(work_surf);
work_surf = work_surf(:,2:size_temp(2));
size_temp = size(work_surf);

psi_coordinates = psi_coordinates';
theta_coordinates = theta_coordinates';

psi = psi_coordinates'; 
th = theta_coordinates';

This 2 snips of code can be modified to remove redundancy

[HN]

Hi @Dyuman Joshi,

The unit is Meridian is rad. Yes, that is right. these lines are redudant and should be removed.

[HN]

@Jeffrey Clark, Right. It woudl be easy to help if the entire working code is uploaded. But I am afraid I cannot do that..