Hi I have this equation with two integration constants

Hi I have this equation with two constants (C1,C2)

and I have these boundary conditions:

u=0 , y=h

u=U , y=0

How I can solve it analytically in Matlab?

I tried to solve it in Maple and I applied the boundary conditions and I got :

and

Thanks

Answer 1

[Walter Roberson]

Are we able to assume G>=0 and that n is a positive integer? Are we able to further assume that you only want real-valued C1 and C2 ?

With substitution, you get

 U = n*G^(1/n)*((-C2*(n+1)*G^(-1/n)/n)^(n/(n+1))-h)^((n+1)/n)/(n+1)+C2

if you solve for C2 then

 C2 = RootOf(-n*G^(1/n)*((-Z*(n+1)*G^(-1/n)/n)^(n/(n+1))-h)^((n+1)/n)-Z*n+U*n-Z+U, Z)

where here RootOf(f(Z),Z) stands for the set of values of Z such that f(Z) = 0 -- the roots of the equation.

But under the assumption that G is non-negative and that n is a positive integer, then the power (n+1)/n can be expanded out. The integer power n+1 can be expanded out using the binomial theorem to find an expanded version of -Z*(n+1)*G^(-1/n)/n and then you can do a change of variables to eliminate the explicit n'th root. And then you can do the same thing again with the -h and power of (n+1)/n .

In other words, for any specific positive integer n, the RootOf corresponds to the root of a polynomial, and in theory the coefficients of that polynomial could be calculated symbolically in terms of the binomial theory applied a couple of times. This all should leave you with roots of a polynomial of potentially high degree.

For example for n = 6 then

C2 = - (RootOf(117649*326592^(1/7)*G^(1/7)*Z^36*h+705894*U*Z^35-50421*326592^(2/7)*G^(2/7)*Z^30*h^2+1764735*U^2*Z^28+12005*326592^(3/7)*G^(3/7)*Z^24*h^3+2352980*U^3*Z^21-1715*326592^(4/7)*G^(4/7)*Z^18*h^4+1764735*U^4*Z^14+147*326592^(5/7)*G^(5/7)*Z^12*h^5-7*326592^(6/7)*G^(6/7)*Z^6*h^6+705894*U^5*Z^7+46656*G*h^7+117649*U^6,Z)^6)^(7/6)

Thus I figure there is an analytic solution for this special case -- but it is not going to be straight forward to express. And "analytic solution" only in the sense of getting down to roots of a polynomial of high degree, which is not to say there is an algebraic solution since the roots will not generally be expressible as radicals.

[Sako Azeez]

Thanks for your Answer actually G>0 and the range of n (0.5 - 1.5) is there any way that I can solve it analytically by iteration ? iterate values for C1 and C2?

[Walter Roberson]

For any rational n in that range, you could perhaps use a similar argument to what I used above, since n/(n+1) of a rational would still be a rational. On the other hand I appear to be having difficulty with, for example, n = 8/11, and it is possible that my argument does not apply.

Is there any way to solve it by iteration? I doubt it with symbolic n. For any fixed numeric n and G and h values, then Yes, it would be possible to solve numerically.

Answer 2

[Abraham Boayue]

Hey Sako, I think Walter has done an excellent job in helping you tackle your equation above. I would like to take a slightly different approach, mainly analytically and a bit of coding. Check the attached file for the former. Follow the written solution before executing the code. An important concept about softwares such as matlab and Mapple is that, they are not good at simplifying expressions like we human do. Mapple was able to give you the correct system of equations that you needed to solve for both C1 and C2, but could not separate the two. Take a look at my solution and let me know what you think.

clear variables
close all
% Conpute the binomial coefficients (For the example considered in the derivation)
  N = 2;
  Sk = [];
  for k = 0:N
      ck = factorial(N)/(factorial(k)*factorial(N-k));
      Sk = [Sk ; k ck];
  end
 %[25 + 10*C1-2 = 0]; Just to illustrate how the root function works.
 % Perhahaps you already know this.
 C1 = roots([0 10 23]);
 disp(C1);

fprintf('----Table of ck -----------\n')
fprintf('    k              ck      \n')
fprintf('===========================\n')
fprintf('%12.5f  %12.5f    \n',Sk');
fprintf('===========================\n');

y = -50:0.1:50;
u = 2*((y-C1).^2-(5-C1)^2);
plot(y,u,'linewidth',2,'color','m');
a = title('u(y)');
set(a,'fontsize',14);
a = xlabel('y[ -50 50]');
set(a,'fontsize',14);
a = ylabel('values'); 
set(a,'fontsize',14);
grid

[Walter Roberson]

This analysis assumes that (n+1)/n is an integer. If we let n=p/q in lowest terms (integers, no common factors) then (n+1)/n is (p/q+1)/(p/q) = (p+q)/q/(p/q) = (p+q)/p . The only time this can be an integer is if q/p is an integer, q=c*p for some integer c. That would make the original p/q into p/(c*p) = 1/c. So n must be 1/c for some integer c. There are only two integers for which 1/c is in the range 1/2 to 3/2, namely c=2 giving n=1/2, and c=1 giving n=1.

Your analysis fails for all other values of n within the required range.

Answer 3

[Abraham Boayue]

Let us note that in my analysis I did not assume anything about n being an integer. Rather, it was the fraction of (1+n)/n that I said should be an integer, which makes perfect sense in that I was able to come up with a simpler solution. Although the information that n should be in the range [ 1/2 3/2] was not available prior to my solution, my method capture two of those values (1/2 and 1 for N = 3 and 2), in fact, for any rational n in the range, we can come up with a corresponding whole number. The proof is easy, but I will only give the result here. I have now revised the solution and made some changes in the code to include only whole numbers within the range of n. I personally believe that it would be very awkward to solve this problem without replacing (1+n)/n by an integer exponent.

clear variables
close all
% Conpute the binomial coefficients 
n0 = 1/2;
n1 = 3/2;
L = 4;
delta = (n1-n0)/L;
n = n0:delta:n1;
Nn = length(n);
N = (1+n)./n;
s = zeros(1,Nn);
M = zeros(1,Nn);
for i = 0: Nn-1
    if ( i == 0)
        s(i+1) = 1;
    else
     s(i+1) = (i/n(1))+(1/delta);
           remainder = rem(s(i+1),2);
            if (remainder == 0)
                s(i+1) = s(i+1)/2;
            end

      end
  end

for i = 1 : Nn

    if (rem(N(i),1)==0 )
           M(i) = N(i);
           s(i) = 1;
      else
         M(i) = s(i)*N(i);
    end

end

n = n';
N = N';
s = s';
M = M';
fprintf('-------------------Table 1 Whole numbers of n --------------\n')
fprintf('       n                  S             N           M = s*N \n')
fprintf('============================================================\n')
fprintf(' %12.5f    %12.5f    %12.5f    %12.5f  \n',  [n'; s'; N'; M']);
fprintf('============================================================\n');

% Compute the binomial coefficients
q =  M(1); % change these in pairs M(2) and N(2)
N1 = N(1);
Sk = [];
  for k = 0:q
      ck = factorial(q)/(factorial(k)*factorial(q-k));
      Sk = [Sk ck];
  end

fprintf('----Table 2 of ck -----------\n')
fprintf('    k              ck      \n')
fprintf('===========================\n')
fprintf('%12.5f  %12.5f    \n',Sk');
fprintf('===========================\n');

h = 5;
G = 4;
U = 4;
Mg = 1/N1*G^(N1-1);
cv = zeros(1, length(Sk));% Coefficient vector for finding C1
cv(1) = Sk(end)-1;
cv(end) = Sk(1)*h^q-U/Mg;
for i = 2:length(Sk)-1

      cv(i) = Sk(i)*h^(i-1);
  end
  %     disp('cv')
  %     disp(cv)
      C1 = roots(cv);
  %     disp('C1')
  %     disp(C1)

% Plot u(y) for real values of C1
y = -50:0.1:50;
u = real(Mg*((y-C1(1)).^N1-(5-C1(1)).^N1));
plot(y,u,'linewidth',2,'color','m');
a = title('u(y)');
set(a,'fontsize',14);
a = xlabel('y[ -50 50]');
set(a,'fontsize',14);
a = ylabel('values'); 
set(a,'fontsize',14);
grid

[Walter Roberson]

"in fact, for any rational n in the range, we can come up with a corresponding whole number."

No. Expand out (1+n)/n . That is 1/n + n/n = 1/n + 1 . The only way this can be an integer is if 1/n is an integer. In the range [1/2, 3/2] the only values that are appropriate are n = 1/2 and n = 1, which lead to N = (1+1/2)/(1/2) = (3/2)/(1/2) = 3, and N = (1+1)/(1) = 2 .

Consider that for N = 4 to be generated then that would require n = 1/3, (1+1/(1/3)) = 1+3 = 4, but that is below the lower bound of the range [1/2, 3/2]

Consider that for N = 1 to be generated then that would require that 1+1/n = 1 would required that 1/n is 0 which would require that n = infinity, which is above the upper bound of the range [1/2, 3/2]

For all other n in the range [1/2, 3/2], (1+n)/n cannot be an integer, so your line that assumes N = (1+n)/n is an integer that can be used for binomial coefficients is wrong except for those two special cases. That ruins the rest of the proof, unfortunately.

Answer 4

[Abraham Boayue]

I get your point and reasoning. The truth is that I am no longer dealing with N directly (N is no longer taken as a whole number), I have scaled N such that it no longer goes above 3 which is the maximum value it can have for our purpose.

If you have the time, I would recommend that you kindly run the code that I posted. Verify that this is true by changing the step size in the range of n. My new variable that will yield all integers within the range is now defined as M = s*N. Where s is a scaled factor that will force N to stay within the range of n.

Pick a value for n, say n = 2/3 which is in range, then N = 1+1/n = 1+3/2 = 5/2, this means that M = s*N = 5, with s = 2. So, as you can see, s is my insurance policy that N will always stay within range.

This example summarizes the basis of my solution, if you can disprove this fact, I might be willing to reconsider my claim.