Symbolic integration: error

Hi to all,

I have been trying to write a piece of code to avoid doing some calculations: the purpose is to calculate the integral of a given function and the integrale of the norm squared of its gradient on a shape that resembles a disk with a point, all of which depend on a parameter r.

For the first test, I did the calculations by hand and everything worked. Then I tried to automate the process of calculating the gradient and its norm squared in two different ways, and in both of them the functions are calcukated but the program fails to calculate the integral.

I've attached both the wrking and non-working tests (sorry for the long lines, but the shape and the integral are not easy to express).

Thank you in advance!

%test1 working
syms x y r

u(x,y,r)=((x-1+r)^2+y^2-r^2)*(y^2-r^2/(1-2*r)*x^2);
du2(x,y,r)=(y^2-r^2/(1-2*r)*x^2)^2*(4*(x-1+r)^2+4*y^2)+((x-1+r)^2+y^2-r^2)^2*(4*(r/sqrt(1-2*r))^4*x^2+4*y^2)-8*r^2/(1-2*r)*x*(x-1+r)*(y^2-r^2/(1-2*r)*x^2)*((x-1+r)^2+y^2-r^2)+8*y^2*(y^2-r^2/(1-2*r)*x^2)*((x-1+r)^2+y^2-r^2);

intu= int( int(u,y,-sqrt(r^2-(x-1+r)^2),sqrt(r^2-(x-1+r)^2)), x, (1-2*r)/(1-r), 1)  + int( int(u,y,-r/sqrt(1-2*r)*x,r/sqrt(1-2*r)*x), x, 0, (1-2*r)/(1-r));
intdu2= int( int(du2,y,-sqrt(r^2-(x-1+r)^2),sqrt(r^2-(x-1+r)^2)), x, (1-2*r)/(1-r), 1)  + int( int(du2,y,-r/sqrt(1-2*r)*x,r/sqrt(1-2*r)*x), x, 0, (1-2*r)/(1-r));

%test 2 not working (first way)
syms x y r

u(x,y,r)=((x-1+r)^2+y^2-r^2)*(y^2-r^2/(1-2*r)*x^2);
du(x,y,r)=gradient(u, [x,y]);
du_body=formula(du);
du2(x,y,r)=du_body(1)^2+du_body(2)^2;

intu= int( int(u,y,-sqrt(r^2-(x-1+r)^2),sqrt(r^2-(x-1+r)^2)), x, (1-2*r)/(1-r), 1)  + int( int(u,y,-r/sqrt(1-2*r)*x,r/sqrt(1-2*r)*x), x, 0, (1-2*r)/(1-r));
intdu2= int( int(du2,y,-sqrt(r^2-(x-1+r)^2),sqrt(r^2-(x-1+r)^2)), x, (1-2*r)/(1-r), 1)  + int( int(du2,y,-r/sqrt(1-2*r)*x,r/sqrt(1-2*r)*x), x, 0, (1-2*r)/(1-r));

%test2 not working (second way)
syms x y r

u(x,y,r)=((x-1+r)^2+y^2-r^2)*(y^2-r^2/(1-2*r)*x^2);
du(x,y,r)=gradient(u, [x,y]);
du2(x,y,r)=(norm(du))^2;

intu= int( int(u,y,-sqrt(r^2-(x-1+r)^2),sqrt(r^2-(x-1+r)^2)), x, (1-2*r)/(1-r), 1)  + int( int(u,y,-r/sqrt(1-2*r)*x,r/sqrt(1-2*r)*x), x, 0, (1-2*r)/(1-r));
intdu2= int( int(du2,y,-sqrt(r^2-(x-1+r)^2),sqrt(r^2-(x-1+r)^2)), x, (1-2*r)/(1-r), 1)  + int( int(du2,y,-r/sqrt(1-2*r)*x,r/sqrt(1-2*r)*x), x, 0, (1-2*r)/(1-r));

[Umar]

To resolve the issues with your integral calculations in MATLAB, it is essential to ensure that the gradient is computed correctly and that the resulting expressions are properly formatted for integration. In your second test, the use of gradient may not yield the expected symbolic output, which can lead to problems when calculating the norm squared.

Answer 1

[Star Strider]

I ran your code in MATLAB Online.

There are some integrands that simply cannot be integrated and give closed-form symbolic results. Yours seem to be among those. One option may be to perhaps do trigonometric substitution (if that is appropriate, since I am not certain what you are doing), and then integrate. If you have actual numeric values instead of symbolic variables, either use integral2 or if you want to stay with the Symbolic Math Toolbos, the vpaintegral function instead of int.

.

[Simone]

Thank you for your answer.

As I understand it, even in the case a closed-form result cannot be given, int should still give produce a result, even though very much implicit, like shown in the help page for int; instead, the program never ends, and when I stop it all the results are there, except for the last integral. Moreover, the function I am differentiating and then integrating is polynomial with respect to those variables, so there shoudl definitely exist a closed-form result...

[Star Strider]

Consider using the simplify function on your expressions, giving the function perhaps 500 iterations for each to see what it can do. (To do the trigonometric substitution, use the subs or rewrite functions.) The int function is clearly having problems with them as they exist.

[Torsten]

If "int" cannot find an antiderivative, it will return the integral expression. But it can take a very long time until MATLAB gives up.

Please explain what you are trying to do - maybe we can find a different solution to your problem.

[Walter Roberson]

For test1, simplify(intu,500) manages to rewrite into an expression roughly 60% the size of the original. However, the int() part of it remains unchanged except for removal of a multiplication by 2.

rewrite(intu,VariousWays) leaves intu unchanged for all of the ways that I tried. (More correctly, I extracted the int() portion of intu and rewrite() that, unsucessfully.)

[Simone]

My goal is mainly this: I have a shape optimization problem, and right now it involves a figure which is a circle centered at (1-r,0) with radius r, united with a sort of triangle with two vertices on the circumference (symmetric wrt the x-axis) and the third in (0,0), with r being a parameter varying in [1-1/sqrt(3),1/2). I have to find a function which is zero at the boundary of this shape, and to begin with I chose the simplest, which is the multiplication of the equations of the circumference and the two straight lines; then I have to calculate the integral of the function on the shape, and the integral of the norm squared of the gradient of the function again on the shape. In the end, I would like to plot the function intu^2/intdu2 with respect to the parameter r.

What I really can't understend is why in the first case, when I calculated the gradient and the norm squared by hand, everything seems to work, but when I tried to automate the process (and I verified that those commands did exactly what I wanted them to do), the script seems to loop on the calculation of the last integral, when in my head it should be the same.

Thank you all or your answers!

[Torsten]

What I really can't understend is why in the first case, when I calculated the gradient and the norm squared by hand, everything seems to work,

But "int" already hangs on your first integral and cannot solve it (and I guess it's the same for all other integrals you try to compute):

syms x y r
u(x,y,r)=((x-1+r)^2+y^2-r^2)*(y^2-r^2/(1-2*r)*x^2);
intu= int( int(u,y,-sqrt(r^2-(x-1+r)^2),sqrt(r^2-(x-1+r)^2)), x, (1-2*r)/(1-r), 1) + int( int(u,y,-r/sqrt(1-2*r)*x,r/sqrt(1-2*r)*x), x, 0, (1-2*r)/(1-r));

As soon as the expression that "int" returns contains an integral expression, "int" wasn't able to get a solution.

I have to find a function which is zero at the boundary of this shape, and to begin with I chose the simplest, which is the multiplication of the equations of the circumference and the two straight lines

The simplest is f = 0. So what other constraints do qualify a function to solve your problem ?

[Simone]

Yes, sorry, I wasn't too specific. I need to find a non-zero function that has zero boundary values and is in H^1 of the shape, so polynomials like those I tried should work.

And yes, I know "int" does not return a closed solution, but still the solution with the integral expression is enough to plot the function intu^2/intdu2; instead, the last integral does not even produce a solution with an integral expression, but just loops.

I hope I clarified, in case let me know and I will try to be more specific.

[Torsten]

the last integral does not even produce a solution with an integral expression, but just loops.

The norm function produces absolute values. Try

syms x y r real

and see if this makes a difference.

Any reason why you don't use "integral" or "integral2" for a predefined vector of numerical values for r ?

[Simone]

Honestly, I hadn't thought about it, but I'd rather not have a discretization of the values I'm looking for...

[Torsten]

I hadn't thought about it, but I'd rather not have a discretization of the values I'm looking for...

If "int" doesn't succeed in computing the antiderivatives, it will use a numerical approach for plotting anyway. And the numerical functions are very much faster than the symbolic ones. So if I were you, I'd give it a try.

Answer 2

[Torsten]

Here is one possible implementation:

syms x y r
u = ((x-1+r)^2+y^2-r^2)*(y^2-r^2/(1-2*r)*x^2);
du = gradient(u, [x,y]);
du2 = du(1)^2+du(2)^2;
u = matlabFunction(u);
du2 = matlabFunction(du2);

r = 0.45;
fun1 = @(r,x)integral(@(y)u(r,x,y),-sqrt(r^2-(x-1+r).^2),sqrt(r^2-(x-1+r).^2));
result1 = integral(@(x)fun1(r,x),(1-2*r)/(1-r),1,'ArrayValued',1);
fun2 = @(r,x)integral(@(y)u(r,x,y),-r/sqrt(1-2*r)*x,r/sqrt(1-2*r)*x);
result2 = integral(@(x)fun2(r,x),0,(1-2*r)/(1-r),'ArrayValued',1);
result = result1 + result2

result = 0.0417

fun1 = @(r,x)integral(@(y)du2(r,x,y),-sqrt(r^2-(x-1+r).^2),sqrt(r^2-(x-1+r).^2));
result1 = integral(@(x)fun1(r,x),(1-2*r)/(1-r),1,'ArrayValued',1);
fun2 = @(r,x)integral(@(y)du2(r,x,y),-r/sqrt(1-2*r)*x,r/sqrt(1-2*r)*x);
result2 = integral(@(x)fun2(r,x),0,(1-2*r)/(1-r),'ArrayValued',1);
result = result1 + result2

result = 0.2009