How can I Plot a primitive function?

I'd like to use ezplot to draw a function. This is my code:

ray=0.11;
m=3.5;
C=2*10^(-12);
DS=140;
CONST=1/((DS^m)*C*pi^(m/2));
syms x
F= @(x) CONST*((x.^(-m/2))./((1.6906*(x/ray).^5-4.6813*(x/ray).^4+5.4253*(x/ray).^3-2.8213*(x/ray).^2+0.8716*(x/ray)+0.4251).^m));
syms x
z=int(F,x,0,0.04);
syms x
ezplot(z,[0,0.04])

I always find an error while I want to plot this function... Why? This is the error :

Error using inlineeval (line 14)
Error in inline expression ==> int(2000./(x.^(7./4).*((2179.*x)./275 - (28213.*x.^2)./121 + (5425300.*x.^3)./1331 - (468130000.*x.^4)./14641 +
(16906000000.*x.^5)./161051 + 4251./10000).^(7./2)), x, 0, 1./25)
 Undefined function 'int' for input arguments of type 'double'.

Error in inline/feval (line 33)
        INLINE_OUT_ = inlineeval(INLINE_INPUTS_, INLINE_OBJ_.inputExpr, INLINE_OBJ_.expr);

Error in ezplotfeval (line 51)
    z = feval(f,x(1));

Error in ezplot>ezplot1 (line 468)
    [y, f, loopflag] = ezplotfeval(f, x);

Error in ezplot (line 144)
                [hp, cax] = ezplot1(cax, f{1}, vars, labels, args{:});

Error in sym/ezplot (line 76)
      h = ezplot(fhandle(f),varargin{:});

Answer 1

[Walter Roberson]

xvals = linspace(0, 0.04, 1000);
y = F(xvals);
z = cumtrapz(y);
plot(xvals, z);

[Nicolas Kennis]

Walter, thanks for your answer. I'm not sure to understand... Can you please write down the all code?

This is what I understand :

 ray   = 0.11;
 m     = 3.5;
 C     = 2*10^(-12);
 DS    = 140;
 CONST = 1/((DS^m)*C*pi^(m/2));
 F     = @(x) CONST*((x.^(-m/2))./((1.6906*(x/ray).^5-4.6813*(x/ray).^4+5.4253*(x/ray).^3- ...
             2.8213*(x/ray).^2+0.8716*(x/ray)+0.4251).^m));
 z     = integral(F,0.0007,0.047);
 xvals = linspace(0.0007, 0.047, 1000);
 y = F(xvals);
 z = cumtrapz(y);
 plot(xvals, z);

But like this, are you sure that I have the points for my primitive function z (and not the points from my function F)?

Many thanks

[Walter Roberson]

 ray   = 0.11;
 m     = 3.5;
 C     = 2*10^(-12);
 DS    = 140;
 CONST = 1/((DS^m)*C*pi^(m/2));
 F     = @(x) CONST*((x.^(-m/2))./((1.6906*(x/ray).^5-4.6813*(x/ray).^4+5.4253*(x/ray).^3- ...
             2.8213*(x/ray).^2+0.8716*(x/ray)+0.4251).^m));
 xvals = linspace(0.0007, 0.047, 1000);
 y = F(xvals);
 z = cumtrapz(y);
 plot(xvals, z);

This is numeric integration.

[Nicolas Kennis]

Walter, Sorry to come back with this topic but I'm not really happy with the results I get with your solution. It doesn't give a good plot because when you change the interval (xvals=linspace(0.0007,0.047,1000)) for another beginning (for example 0.001), the plot is also changing. It was the reason I need the primitive function (antiderivative) and not a numerical answer. Is it possible to have something like :

 ray   = 0.11;
m     = 3.5;
C     = 2*10^(-12);
DS    = 140;
CONST = 1/((DS^m)*C*pi^(m/2));
F     = @(x) CONST*((x.^(-m/2))./((1.6906*(x/ray).^5-4.6813*(x/ray).^4+5.4253*(x/ray).^3- ...
             2.8213*(x/ray).^2+0.8716*(x/ray)+0.4251).^m));
z=int(F,x);
plot(z,...)
....

? I insist with the 'int(F,x)' because it's really important for me to get the exact plot and not an approximation.

Thanks in advance

Nicolas

[Walter Roberson]

You cannot get an exact plot. No display can position points with infinite precision. You cannot position a point at x = 1, y = Pi exactly. You can position a point at x = 1, y = some approximation of Pi.

There is no closed form formula for your function, Numeric integration is your only choice. Numeric integration is always approximate. The better numeric integration routines are adaptive, but they all come down to trapezoidal rule or similar applied with respect to the y values at x positions they sampled at.

If you want more precision, sample more finely. Sample at irregular intervals if you want. logspace() instead of linspace().

Your formula for F has only x as its free variable, and in your original post you defined z=int(F,x,0,0.04); which is definite integration of a function with respect to a single variable . The result is going to be a specific numeric value. But then you want to plot that with respect to x.

Your more recent posting uses indefinite integration for z, and asks to plot the result, but it does not define either a left or right boundary, so there cannot be any solution.

Your formula for F approaches infinity quickly as x approaches 0 from the right. The exact value of the integral is going to depend strongly on where you start integrating -- and in turn that means that the plot is going to depend a lot on where you start integrating. It is not the case that there is "thin" infinity were the values might be large but very narrow in such a way that the integral stays finite. Each additional decimal place closer to 0 gives an F value about 56 times larger (10^(7/4) is the ratio)

[Nicolas Kennis]

Walter, here you'll find the code I'm using:

 ray   = 0.11;
m     = 3.5;
C     = 2*10^(-12);
DS    = 140;
CONST = 1/((DS^m)*C*pi^(m/2));
F     = @(x) CONST*((x.^(-m/2))./((1.6906*(x/ray).^5-4.6813*(x/ray).^4+5.4253*(x/ray).^3- ...
             2.8213*(x/ray).^2+0.8716*(x/ray)+0.4251).^m));
F2     = @(x) CONST*((x.^(-m/2))./((2.3652*(x/ray).^5-7.1218*(x/ray).^4+8.5106*(x/ray).^3- ...
             4.6314*(x/ray).^2+1.5432*(x/ray)-0.0088).^m));
xvals = linspace(0.0007, 0.003, 1000);
xvals2 = linspace(0.0007, 0.003, 1000);
y = F(xvals);
y2 = F2(xvals2);
z = cumtrapz(y);
z2 = cumtrapz(y2);
plot(xvals, z,'color','r');hold on;
plot(xvals2, z2,'color','b');
view(-90,90);
set(gca,'ydir','reverse')
xlabel('Profondeur de fissure (m)');
ylabel('Nombre de cycles (N)');
set(gca,'YScale','log');
grid on

It gives me a result. When I tried to put logspace() instead of linspace(), I've have a very big difference. Is it normal? It seems to be good with linspace() but not with logspace()...

[Walter Roberson]

When you use logspace remember that the arguments should be the power of 10... eg, logspace(-1,2) instead of logspace(.1,100)

[Walter Roberson]

You forgot to pass the x values to cumtrapz()

[Nicolas Kennis]

Yes, I saw the solution for logspace() for the power of 10. But I don't understand what you mean with cumtrapz(), I followed your example to write my code...? Could you please write down the good code? This is what I think :

 ray   = 0.11;
m     = 3.5;
C     = 2*10^(-12);
DS    = 140;
CONST = 1/((DS^m)*C*pi^(m/2));
F     = @(x) CONST*((x.^(-m/2))./((1.6906*(x/ray).^5-4.6813*(x/ray).^4+5.4253*(x/ray).^3- ...
             2.8213*(x/ray).^2+0.8716*(x/ray)+0.4251).^m));
F2     = @(x) CONST*((x.^(-m/2))./((2.3652*(x/ray).^5-7.1218*(x/ray).^4+8.5106*(x/ray).^3- ...
             4.6314*(x/ray).^2+1.5432*(x/ray)-0.0088).^m));
xvals = logspace(-3.5, -2.5);
xvals2 = logspace(-3.5, -2.5);
y = F(xvals);
y2 = F2(xvals2);
z = cumtrapz(y);
z2 = cumtrapz(y2);
plot(xvals, z,'color','r');hold on;
plot(xvals2, z2,'color','b');
view(-90,90);
set(gca,'ydir','reverse')
xlabel('Profondeur de fissure (m)');
ylabel('Nombre de cycles (N)');
set(gca,'YScale','log');
grid on

[Walter Roberson]

z = cumtrapz(xvals, y);
z2 = cumtrapz(xvals2, y2);

[Nicolas Kennis]

Finally, I have this:

 ray   = 0.11;
m     = 3.5;
C     = 2*10^(-12);
DS    = 140;
CONST = 1/((DS^m)*C*pi^(m/2));
F     = @(x) CONST*((x.^(-m/2))./((1.6906*(x/ray).^5-4.6813*(x/ray).^4+5.4253*(x/ray).^3- ...
             2.8213*(x/ray).^2+0.8716*(x/ray)+0.4251).^m));
F2     = @(x) CONST*((x.^(-m/2))./((2.3652*(x/ray).^5-7.1218*(x/ray).^4+8.5106*(x/ray).^3- ...
             4.6314*(x/ray).^2+1.5432*(x/ray)-0.0088).^m));
xvals = logspace(-3.5, -2.5);
xvals2 = logspace(-3.5, -2.5);
y = F(xvals);
y2 = F2(xvals2);
z = cumtrapz(xvals,y);
z2 = cumtrapz(xvals2,y2);
plot(xvals, z,'color','r');hold on;
plot(xvals2, z2,'color','b');
view(-90,90);
set(gca,'ydir','reverse')
xlabel('Profondeur de fissure (m)');
ylabel('Nombre de cycles (N)');
set(gca,'YScale','log');
grid on;

Is it correct?

[Nicolas Kennis]

Hello, I do evolved in my function and I've a new code. I've an error but I don't know why? Can you help me? My goal is tho find this integral! I have this error : Error using * , Inner matrix dimensions must agree.

 ray   = 0.11;
m     = 3.5;
C     = 2*10^(-12);
DS    = 140;
CONST = 1/(C*pi^(m/2));
F     = @(x) CONST*((x.^(-m/2))./(((452092178*x.^3-10493093*x.^2+78093.85*x-42.08).^m).*((1.6906*(x/ray).^5-4.6813*(x/ray).^4+5.4253*(x/ray).^3-2.8213*(x/ray).^2+0.8716*(x/ray)+0.4251).^m)));
W=integral(F,0.002,0.01);
disp(['intégrale au point A = : ' num2str(W)]);

[Nicolas Kennis]

OK I found the problem... ;-)

Answer 2

[Mischa Kim]

Nicolas, I would assume that there is no closed-form, symbolic solution for your integral. So instead use numeric integration and the integral function

 ray   = 0.11;
 m     = 3.5;
 C     = 2*10^(-12);
 DS    = 140;
 CONST = 1/((DS^m)*C*pi^(m/2));
 F     = @(x) CONST*((x.^(-m/2))./((1.6906*(x/ray).^5-4.6813*(x/ray).^4+5.4253*(x/ray).^3- ...
             2.8213*(x/ray).^2+0.8716*(x/ray)+0.4251).^m));
 z     = integral(F,0.01,0.04)

Also, because of the first term in F and the lower limit of x the integral will be infinite. I changed it to 0.01 for a finite result.

[Nicolas Kennis]

Kim, Thanks for your answer! I really appreciate.

So, it means that there is no solution to plot the primitive of F (from 0.0007 to 0.05 for example)?

Thanks in advance