Symbolic integration in matlab

I have: function problemOne

%declare symbolic Variables
N=6.023e23;
k=1.3806503e-23; %units (m^2)kg/(s^2)K
e=k*213.96; %units of k * 1/K
r0 = 426; %units of pm

syms r T;

%set y to be function inside integral in B
y=exp(((-4*e)*(((r0/r)^12)-((r0/r)^6))/k*T)-1)*(4*pi*(r^2));

%define B in terms of y
B= (-1/2)*N*int(y,r,0,Inf);

%solve for B=0;

solve(B,T)

end

This outputs:

Warning: Explicit solution could not be found. 
> In solve at 81
  In problemOne at 21

ans =

[ empty sym ]

I think the problem is in the integral, since

int(y,r,0,Inf)

yields

Inf - limit(4*pi*r^3*exp(-1)*exp(4005180207889834777988627621265408/783018253290173*T/r^6)*exp(-23937547659622311621519348000222281120506532855808/783018253290173*T/r^12), r = 0, Right)

What do I do? I thought I had a basic knowledge of Matlab (not an expert by any means), but this is confusing me to no end.

Thanks, Nkk

Answer 1

[William]

When you're trying to run the whole thing at once debugging can be next to impossible. My guess is you're getting infinity somewhere in the problem or you have a syntax error within the code.

Try this.

Break it into pieces. have the program print out all the pieces and check the answers to see if they mathematically make sense.

I cannot see what you have specifically written on line 21 which makes it difficult to see what is going on.

Also see if you can plot the integral before you attempt to integrate it. If it is going to infinity than integrating it will lead to infinity.

try the

ezplot() function.

Answer 2

[Walter Roberson]

y goes to infinity as r goes to infinity. The exp() part has constants divided by r^6 and r^12, so those divisions are going to go to 0 quickly (regardless of T), leading to exp(0). The r^2 term at the end of y therefor dominates, and r^2 strictly increases as r increases so the behavior of y going to infinity is not "removable".

The value at the other end, near r = 0, depends a great deal upon the value of T, especially when r is near 426 and T is near 0.

Answer 3

[nathan khosla]

I know...I think that is my problem. However, I have no idea how to solve it. In the end, I really only want to know when B=0. Is there any easier way to do this? Do I need symbolics? I have a friend who is a bit better at Matlab than I am (read:he is very good) and he said that I may not need them (it was over the phone, so he never saw the code, thus the uncertainty). So...do I need them?

Thanks, Nkk

[Walter Roberson]

As B is a constant multiplied by the integral, B will be 0 if and only if the integral is 0. As indicated above, though, the integral goes to infinity as r increases, and as that is not over an infinitely small region, the integral will always be infinite for every value of T, even for T = 0 or T negative.

As at least some points are going to come out positive, in order to get 0 as the integral you would have to have some points come out negative. exp() is, however, never negative unless T is complex. (And even if T is complex, the division by r^6 or r^12 sends the exponential quickly.)

The only remaining possibility is if T was complex and the integral in the region near 0 was so large as to cancel out the remaining infinite tail. Is a complex T plausible in the situation?