How to convert a symbolic output to its numerical form ( for vector integration) ?

I am stuck in a resilient difficulty due to my inexperience with vector computation.

The code is all about calculating outward vector integrals on a parameterized surface.

The parametric form I have used here is a superellipsoid, which has a general form

(x^eps1)/a1+(y^eps2)/a2+(z^eps3)/a3=1.

With suitable parameterization, I have formed a matrix (x, y, z) which returns the geometry for superellipsoidal shapes, as,

function [x,y,z]=superellipse(epsilon1,epsilon2,epsilon3,a1,a2,a3)
  n=100;
  etamax=pi/2;
  etamin=-pi/2;
  wmax=pi;
  wmin=-pi;
  deta=(etamax-etamin)/n;
  dw=(wmax-wmin)/n;
  [i,j] = meshgrid(1:n+1,1:n+1);
  eta = etamin + (i-1) * deta;
  w   = wmin + (j-1) * dw;
  x = a1.*sign(cos(eta)).*abs(cos(eta)).^epsilon1.*sign(cos(w)).*abs(cos(w)).^epsilon1;
  y = a2.*sign(cos(eta)).*abs(cos(eta)).^epsilon2.* sign(sin(w)).*abs(sin(w)).^epsilon2;
  z = a3.*sign(sin(eta)).*abs(sin(eta)).^epsilon3;
  surf(x,y,z);
  surf2stl('cube3D1.stl',x,y,z,'binary');
  end

now for the second part, we need to compute the pressure tensor on the surface points. I follow the procedure of integral mentioned here. Please follow the link.

syms eta w real
ellip=[(a1.*sign(cos(eta)).*abs(cos(eta)).^epsilon1.*sign(cos(w)).*abs(cos(w)).^epsilon1),(a2.*sign(cos(eta)).*abs(cos(eta)).^epsilon2.*sign(sin(w)).*abs(sin(w)).^epsilon2),(a3.*sign(sin(eta)).*abs(sin(eta)).^epsilon3)];
F=[(a1.*sin(eta).*cos(w)),(a2.*sin(eta).*sin(w)),(a3.*cos(eta))];
nds=simplify(cross(diff(ellip,eta),diff(ellip,w)));hold on
Fpar=subs(F,[(a1.*sin(eta).*cos(w)),(a2.*sin(eta).*sin(w)),(a3.*cos(eta))],ellip);
Fj=@(Fpar,nds)Fpar.*transpose(nds);
flux=(symint2(Fj,eta,-pi,pi,w,-(pi/2),(pi/2)));

I have solved the integral analytically with Matlab syms function. But, after computing the integral,

the values have been returned in 'char' form. I need to have the values in 'double'.

I have used a function symint2, which is attached herewith.

Please take a look at the code.

Hirak

Answer 1

[madhan ravi]

Simply use double() or vpa() to convert symbolic answer to numerical also have a look at matlabFunction() to convert symbolic operations to numeric ones thereby you can use numerical integration (integral()) too.

[Hirak]

Please elaborate.

I have tried to put vpa() on nds. doesnt result any fruitful solution.

Thanks for mentioning matlabFunction() .I will definitely try to resolve with it.

[madhan ravi]

Anytime :) , the code isn't organised so having trouble with understanding your code , please try the changes with matlabFunction() and get back if any confusions

[Hirak]

I have tried but failed to use matlabFunction to define Fj in the code. The vpaintegral also does not wok. please help

%%%%%%%%%%%%%%% defigning superellipsoid %%%%%%%%%%%%%%%%%
 %declare constants%
 a1=25; a2=25; a3=25; epsilon1=1; epsilon2=1; epsilon3=1;n=100;

 %declare variables%
 etamax=pi/2; etamin=-pi/2;
 wmax=pi; wmin=-pi;
 deta=(etamax-etamin)/n;
 dw=(wmax-wmin)/n;
 [i,j] = meshgrid(1:n+1,1:n+1);
 eta = etamin + (i-1) * deta; 
 w   = wmin + (j-1) * dw;
 x = a1.*sign(cos(eta)).*abs(cos(eta)).^epsilon1.*sign(cos(w)).*abs(cos(w)).^epsilon1;
 y = a2.*sign(cos(eta)).*abs(cos(eta)).^epsilon2.* sign(sin(w)).*abs(sin(w)).^epsilon2;
 z = a3.*sign(sin(eta)).*abs(sin(eta)).^epsilon3;
 shape=surf(x,y,z);
 ellipse=[x,y,z];
 
 %ellipse defines the set of points on the surface of ellipsoid
 
 %%%%%%%%%%%%%%%%%%%% defining force %%%%%%%%%%%%%%%%%%%%%%%%
 
 x1=a1.*sin(eta).*cos(w);
 y1=a2.*sin(eta).*sin(w);
 z1=a3.*cos(eta);
 F=[(a1.*sin(eta).*cos(w)),(a2.*sin(eta).*sin(w)),(a3.*cos(eta))];        

 %%%%%%%%% parameterizing values of F on ellipse to satisfy one another (using symbolic variables)%%%%%%%%%
 syms eta w
 ellip_par=[(a1.*sign(cos(eta)).*abs(cos(eta)).^epsilon1.*sign(cos(w)).*abs(cos(w)).^epsilon1),(a2.*sign(cos(eta)).*abs(cos(eta)).^epsilon2.*sign(sin(w)).*abs(sin(w)).^epsilon2),(a3.*sign(sin(eta)).*abs(sin(eta)).^epsilon3)];       
 F_par=subs(F,[(a1.*sin(eta).*cos(w)),(a2.*sin(eta).*sin(w)),(a3.*cos(eta))],ellip_par);
 
 %F_par defines the set of points on to which the force is applied to the
 %surface of ellipsoid, satisfying each elements of predefined matrices

 %%% now, we have to evaluate the surface vector integral poiting outside
 %%% to the surface to evaluate total pressure on the surface. This has to
 %%% in coupled of steps. First, we need to define the outward vector matrix.
 
 %%%%%%%%%%%%%%%%%%%% Defining the outward vector matrix %%%%%%%%%%%%%%%%%%%
 nds=cross(jacobian(ellip_par,eta),jacobian(ellip_par,w));
 Fj=@(Fpar,nds)Fpar.*transpose(nds);
 %Fparam=matlabFunction(Fj);
 flux=vpaintegral(Fj(eta,w),eta,-pi,pi,w,-(pi/2),(pi/2));

[madhan ravi]

lot of mess:

 Fj=matlabFunction(F_par.*transpose(nds));
 %Fparam=matlabFunction(Fj);
 flux=integral(Fj,-pi,pi,-(pi/2),(pi/2));

• size(F_par) is 101 X 303
• size(transpose(nds)) is 1 X 3

Now you have to decide how you want to multiply F_par with transpose(nds) (it's simply nds.' ).

[Hirak]

Dear Sir,

Thank you for your critical reply.

I have understood my problem and in order to solve that, I have divided the code in two parts. first one is to compute the outgoing flux

function[flux]= superquad(a1,a2,a3, epsilon1,epsilon2,epsilon3,eta,w)
 syms eta w x y z
 elp_x = a1.*sign(cos(eta)).*abs(cos(eta)).^epsilon1.*sign(cos(w)).*abs(cos(w)).^epsilon1;
 elp_y = a2.*sign(cos(eta)).*abs(cos(eta)).^epsilon2.* sign(sin(w)).*abs(sin(w)).^epsilon2;
 elp_z = a3.*sign(sin(eta)).*abs(sin(eta)).^epsilon3;
 ellipsoid=[elp_x,elp_y,elp_z];
 F=[x,y,z];
 nds=cross(diff(ellipsoid,eta),diff(ellipsoid,w));
 F_par=subs(F,[x,y,z],ellipsoid);
 realdot=@(u,v)u*transpose(v);
 flux=symint2(realdot(F_par,nds),w,-pi/2,pi/2,eta,-pi,pi);
end

Which returns me the functional form very well. For the next part I would like to calculate the output through a for loop, something like that:

%%%%%%%%%%%%%%% defigning superellipsoid %%%%%%%%%%%%%%%%%
%declare constants%
 a1=25; a2=25; a3=25; epsilon1=1; epsilon2=1; epsilon3=1;n=100;
 
%declare constants%
 thetamax=pi/2; thetamin=-pi/2;
 phimax=pi; phimin=-pi;
 dtheta=(thetamax-thetamin)/n;
 dphi=(phimax-phimin)/n;
 [i,j] = meshgrid(1:n+1,1:n+1);
 theta = thetamin + (i-1) * dtheta; 
 phi= phimin + (j-1) * dphi;
 
%invoking superquadrics%
for theta=i
     i=i+1;
     for phi=j
         j=j+1;
 superquad(a1,a2,a3,epsilon1,epsilon2,epsilon3,theta,phi);
     end
end

In this case, solution did not come after a prolonged run.

Please help me.

[madhan ravi]

You can unaccept my answer so that someone else can help , I am not able to figure it out.

[Hirak]

ok.....

Answer 2

[Walter Roberson]

 flux=vpaintegral(Fj(eta,w),eta,-pi,pi,w,-(pi/2),(pi/2));

vpaintegral() does not have a syntax for 2D integrals. You need

  flux=vpaintegral(vpaintegral(Fj(eta,w),eta,-pi,pi),w,-(pi/2),(pi/2));
  

This will give you an answer that is 0 to within round-off.

This is because...

  int(Fj(eta,w),eta,-pi,pi)

is 0 becaus Fj(eta,w) is simply eta*w

You need to look more closely at your lines

 nds=cross(jacobian(ellip_par,eta),jacobian(ellip_par,w));
 Fj=@(Fpar,nds)Fpar.*transpose(nds);
 %Fparam=matlabFunction(Fj);
 flux=vpaintegral(Fj(eta,w),eta,-pi,pi,w,-(pi/2),(pi/2));

The first of those lines defines a variable named nds . The second line defines an anonymous function with dummy parameter name nds that does a multiplication involving the dummy variable nds . With nds being a dummy variable name, the nds after the @(Fpar,nds) has nothing to do with the variable nds that was define in the previous line.

Now, the third line invokes Fj(eta,w) with the symbolic variables eta and w. Those get passed into Fj(), where they become locally known as Fpar and nds . transpose(nds) is then transpose(w) which is just w. Fpar is the passed in eta. So Fj(eta,w) is going to be just eta .* w .

I have modified your code to get closer, but I am not convinced that this is what you are looking for.

Answer 3

[Hirak]

Dear Sir,

I have made some changes to avoid the unattended complexities arising due to the complex nature of functional forms. As it was mentioned in earlier comments, I have defined the vector integral separately through a matlab function, which returns me the symbolic result of integration. The code is mentioned below.

function[flux]= superquad2(a1,a2,a3, epsilon1,epsilon2,epsilon3)
 syms eta w x y z 
 wmin=-pi/2;
 wmax=pi/2;
 etamin=-pi;
 etamax=pi;
 elp_x = a1.*abs(cos(eta)).^epsilon1.*abs(cos(w)).^epsilon1.*sign(cos(eta)).*sign(cos(w));
 elp_y = a2.*abs(cos(eta)).^epsilon2.*abs(sin(w)).^epsilon2.*sign(cos(eta)).*sign(sin(w));
 elp_z = a3.*abs(sin(eta)).^epsilon3.*sign(sin(eta));
 ellipsoid=[elp_x,elp_y,elp_z];
 ndseta=diff(ellipsoid,eta);
 ndsw=diff(ellipsoid,w);
 nds=cross(ndseta,ndsw);
 F=[(a1.*abs(cos(eta)).*abs(cos(w))),(a2.*abs(cos(eta)).*abs(sin(w))),(a3.*abs(sin(eta)))];
 F_par=subs(F,[(a1.*sign(cos(eta)).*abs(cos(eta)).*sign(cos(w)).*abs(cos(w))),(a2.*sign(cos(eta)).*abs(cos(eta)).*sign(sin(w)).*abs(sin(w))),(a3.*sign(sin(eta)).*abs(sin(eta)))],ellipsoid);
 realdot=@(u,v)u*transpose(v);
 fflux=realdot(F_par,nds);
 flux=symint2(fflux,w,wmin,wmax,eta,etamin,etamax);
end

Now I am stuck in the formulation of the second step. I want to compute the numerical values of the symbolic answer for the set of grid points which define the shape of superquadrics, which is defined by the two variables viz

eta

and

w

This part generates the list of gridpoints

%%%%%%%%%%%%%%% defigning superellipsoid %%%%%%%%%%%%%%%%%
%declare constants%
 a1=25; a2=25; a3=25; epsilon1=1; epsilon2=1; epsilon3=1;n=100;
 
%declare constants%
 etamax=pi/2; etamin=-pi/2;
 wmax=pi; wmin=-pi;
 deta=(etamax-etamin)/n;
 dw=(wmax-wmin)/n;
 [i,j] = meshgrid(1:n+1,1:n+1);
 theta = thetamin + (i-1) * dtheta; 
 phi= phimin + (j-1) * dphi;

I tried to feed it all in a for loop like

for theta=i
     i=i+1;
     for phi=j
         j=j+1;
 superquad2(a1,a2,a3,epsilon1,epsilon2,epsilon3,theta,phi);
     end
end

But it continues running without giving a solution. After running an hour I called off the job. Probably its not working. What should I do now? How can I couple the both ends together? Please advise.

Wishing a very prosperous new year,

Hirak

[Walter Roberson]

We do not have a value for thetamin or dtheta or phimin so we cannot execute

 theta = thetamin + (i-1) * dtheta;
 phi= phimin + (j-1) * dphi;

On the other hand your loops after that go

 for theta=i
     i=i+1;
     for phi=j
         j=j+1;
 superquad2(a1,a2,a3,epsilon1,epsilon2,epsilon3,theta,phi);
     end
end

which overwrite theta and phi with entries from the results of the meshgrid() . After the meshgrid, i and j will be 2D arrays of results. When you use for theta=i and i is 2D then the result is defined as assigning theta to be successive columns of i, so the first iteration through theta would be the entire column vector i(:,1), the second iteration it would be the entire column vector i(:,2) and so on. It is unlilkely that is the result you want. Also, you are not assigning the results of the superquad2 call to anything.

[Hirak]

I got your points, and revised the code as;

%%%%%%%%%%%%%%% defigning superellipsoid %%%%%%%%%%%%%%%%%
 a1=25; a2=25; a3=25; epsilon1=1; epsilon2=1; epsilon3=1; n=10;
 
%declare constants%
 etamax=pi/2; etamin=-pi/2;
 wmax=pi; wmin=-pi;
 deta=(etamax-etamin)/(n-1);
 dw=(wmax-wmin)/(n-1);
 for n=1:n
     theta(n)=etamin+(n-1)*deta;
     phi(n)=wmin+(n-1)*dw;
 end 
%invoking superquadrics%
fluxx=superquad2(a1,a2,a3,epsilon1,epsilon2,epsilon3);
for eta= theta(n)
    for w=phi(n)
    flux(eta,w)=vpaintegral(fluxx,eta,etamin,etamax,w,wmin,wmax);
    end
end

Now, the problem appears in running the nested loop invoking the symbolic function. I am failing to substitute the values of eta and w from the array theta and phi, respectively, inside the nested loop defined here.