solve equation with integral with the unknown inside the integral in matlab

please, I am looking for a matlab code which makes it possible to find for what rate rate 'taux1' this equation is verified knowing that Bi and R1 and R2 and n are common

f1(taux1) is:

[ZOUBAIR DAOUMA]

epsilon=-1

[darova]

What about fsolve?

[ZOUBAIR DAOUMA]

the problem is not simple to do just fsolve because rate1 is unknown

[Walter Roberson]

Which variables have known numeric value?

You said before that "Bi and R1 and R2 and n are common" -- does that mean they are all numeric?

Is rate1 the same as R1 ? Are you indicating that it is unknown, and you want to know the solution for as a formula in terms of ?

[ZOUBAIR DAOUMA]

no R1 and n and R2 and B 'are known for example we can give a R1 = 7.40 and R2 = 8.40 and n = 2/3 and B' = 100 and we want for which value of 'taux1' the equation is verified you should also know that the 'taux1' is negative but on not know this values

[Walter Roberson]

Is there a value for epsilon ?

[Walter Roberson]

By the way: given particular R1, R2, B' values and n=2/3, Maple is able to find a closed form integral even without knowing the value of epsilon. However, when I assume epsilon > 0, the equation appears that it would not have a negative tau_1 solution.

[ZOUBAIR DAOUMA]

epsilon= - 1 , please if you can write me the programm matlab if you want

[Walter Roberson]

I assume you need to solve it for multiple parameters. Which values can we assume will not be changed, and which values will be changed?

[Walter Roberson]

I thought I was getting somewhere, but with the values of the constants you posted, there is no real-valued solution for tau_1 . The integral is complex-valued except between approximately -128.8531775 and -100, and everywhere in the range, the value of 1+R1*f1(tau_1)*epsilon is strictly negative (below about -54)

[ZOUBAIR DAOUMA]

for the values ​​that change it is n and R1 and R2 and B 'for epsilon it is always -1 because my goal is to plot this equation (this is the basic speed) if how do you feel please please because i have no idea

g1(r) :

that's why you have to know first 'rate1' so that I can trace the velocity

[ZOUBAIR DAOUMA]

for example you can test n=1; R1=7.40; R2=8.40; Bi=7 and i want to plot Vb so first you have to determine 'taux1' , after I don't know how to plot Vb with the equation I published

[Walter Roberson]

Unfortunately using different n values makes this more difficult as it can change the integral form a lot.

[ZOUBAIR DAOUMA]

I have forget to say that 0 = <n <= 1

[ZOUBAIR DAOUMA]

I know that the integral is not resolved for all value of n but there will exist the cases where it is possible to obtain 'rate1' and to plot Vb

[Walter Roberson]

I find that for any given R1, if n is constrained to be a positive rational and Bprime is constrained to be positive, that I seem to be able to get a closed form solution (in terms of hypergeom function). I am probing further to see if I can get a solution for general positive R1: it has been reluctant to compute one so-far.

[ZOUBAIR DAOUMA]

R1 and R2 and B' its always positif also 0=<n<=1

[Walter Roberson]

Interestingly, the general solution fails for n = 1/q for positive integer q, including for the case n = 1 (which has a solution.)

I also seem to be finding that the overall equation has no solution for at least some ranges of values. There is a relationship that appears to be tau_1 <= - (R1^2*R2^2*Bprime) that must be true for the integral to be real-valued, and for at least some n that is also the minimum of the integral (decreasing tau_1 further increases the integral), but for at least some combinations, the integral at that point is such that 1+epsilon*R1*f1(tau_1) cannot equal 0

[Walter Roberson]

Like 7.4 to 7.5, Bprime = 98, n=2/5 finds a closed form integral in terms of tau_1 in a relatively reasonable time, but 7.4 to 8.4, Bprime = 100, n=2/5 is taking a long time. At the moment I am not sure if the integral can be real-valued with Bprime above 98.3246 for n=2/5

[ZOUBAIR DAOUMA]

I agree with you the integral is not resolable for all the values ​​but I will be so that it values ​​possible to have value of 'rate1' if you please is what you can pass me the matlab coode to check that, if you will

[David Goodmanson]

Hi Zoubair,

Going back to the original question, and denoting eps*taux = etaux (a positive quantity), it looks like you want to find a taux such that

f1 = 1/R1

where f1 is an integral that depends on taux,R1,R2,Bi,n and the last four of these are provided. One thing you can do, at least for the integral calculation, is reduce the number of parameters from four to three. With r = R1*u, dr = R1*du you arrive at

f1 = Int{1,R21} (etaux/u^(n+2) - Bi/u^n)^(1/n) du   where R21 = R2/R1.

There are still four parameters overall since the first equation depends on R1 and not R21.

The integrand is real only when etaux >= etauxmin = R21^2*Bi. For your numbers below

R1 = 7.4;
R2 = 8.4;
Bi = 100;
n = 2/3;
R21 = R2/R1;
etaux_min = R21^2*Bi
etaux = etaux_min     % for function input purposes

fun = @(u) (etaux./u.^(n+2) -Bi./u.^n).^(1/n);
f1_min = integral(fun,1,R21)

f1_min = 7.4355

It's very important that as a function of etaux, f1 is monotonically increasing. You are trying for

f1 = 1/R1 = .1351.

But for the minimum etaux, f1 is already much greater than this, so a solution in the reals is not possible. There are many ways to adjust parqameters to get there but if, for example, you change Bi to 6 and rerun this, then etaux_min is much smaller,

f1_min = .1093 

and there will be a solution. So there is a yes/no test for proposed parameters, although it's a bit inconvenient.

Suppose you want f1 = A, and you have an etaux_min that produces an f1_min that is smaller than A. If you want to solve for etaux by using, say, the fzero function, then etaux_min is obviously a lower search limit and, although it might possibly be on the large side,

etaux_upper = etau_min*(A/f1_min)^n

is guaranteed to give a value of f1 that is greater than A. So etaux_upper can be the upper search limit.

[David Goodmanson]

Hello Zoubair,

to get the velocity you can use, for example

r = R_lower:.001:R2;     % use lots of points;  R_lower might be R1
r = r - r(end) + R2;
G = (etaux*R1^2./r.^(n+2) - Bi./r.^n).^(1/n);
g1 = cumtrapz(r,G);
g1 = g1(end) - g1;             % g1 = 0 at upper limit
V = r.*g1;

I used a convenient spacing for the r array and adusted it to hit R2 exactly, but possibly not R_lower. There are other choices, such as linspace(R_lower,R2,npoints), which might or might not have convenient spacing.

I'm curious about the units of the variables. Have all the variables been normalized from the beginning to be dimensionless?

Ignoring the presence of epsilon for the moment, from the f1 integral, taux and Bi have the same units. f1 and g1 definitely have the same units. But if V is velocity, then V_b (r) = r*g1 seems inconsistent with 1 =R1*f1 (unless all the variables were normalized to be dimensionless at the beginning). If dimensions were restored, f1 and g1 would have to have units of inverse time, and taux would have funny but possible units in terms of length and time.

[ZOUBAIR DAOUMA]

really thank you very much for your help i will see is what the plot program of Vb it works, for all the variables are dimensionless