Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Integrating over a constant function

I am trying to integrate over a constant function in MATLAB 2017a, but I am stuck. First of all when I integrate using the following script, I get the right output. So the script works for a x0 which depends on t.

function E=sol(n,k)
x0 = @(t)  t^(2);
j = 0;
E = zeros(n,1);
 while j < n+1 ;
    K = matlabFunction(subs(po(j,k))) ; 
    eval(sprintf('x%d = integral(K,0,1);',j+1)) ;
    E(j+1,1) = subs(sprintf('x%d',j+1))
    j = j+1;
 end
end

Where function po(j,k) is as follows,

function A_j = po(j,k)          % Adomian polynomials
 if j >0
  x =  sym('x',[1 j]);
    syms p;                            % Assinging a symbolic variable for p
    syms x0;
    S = x0+ sum(p.^(1:j) .* x) ;     % Sum of p*x up to order j
    Q =f(S,k);                        % Taking the k-th power of S, i.e. 
    A_nc = diff(Q,p,j)/factorial(j);  % Taking the j-th order derivative 
    A_j = subs(A_nc,p,0) ;            % Filling in p=0
 else
    syms x0;
    S = x0; 
    A_j =f(S,k);                      % Taking the k-th power of S,
  end
 end 

And where f(x,k) is,

function F = f(x,k) % Nonlinear function of k power
 F = x^k ;
end

Now when I cal sol(n,k) it does work. But when I try to change my x0 function in sol(n,k) in a constant function like,

function E=solcon(n,k)
x0 = @(t) 2.*ones(size(t));
j = 0;
E = zeros(n,1);
 while j < n+1 ;
    K = matlabFunction(subs(po(j,k))) ; 
    eval(sprintf('x%d = integral(K,0,1);',j+1)) ;
    E(j+1,1) = subs(sprintf('x%d',j+1))
    j = j+1;
 end
end

It does not work, As you can see I added *ones(size(t)); just to make it a function of t. But unfortunately it still doesn't work when I call,

K = matlabFunction(subs(po(j,k))) ;

I get,

@()4.0

And then I get an error when I call,

eval(sprintf('x%d = integral(K,0,1);',j+1)) 

Could anyone help me out trying to integrate over a constant?

The error I get when I call solcon(10,2) is

Error using symengine>@()4.0
Too many input arguments.

Error in integralCalc/iterateScalarValued (line 314)
            fx = FUN(t);

Error in integralCalc/vadapt (line 132)
        [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
    [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in solcon1 (line 7)
eval(sprintf('x%d = integral(K,0,1);',j+1)) ;

EDIT 2 I used the following script,

function E=solcon(n,k)
x0 = @(t) 2.*ones(size(t));
j = 0;
E = zeros(n,1);
 while j < n+1 ;
    K = matlabFunction(subs(po(j,k))) ; 
    fstr= func2str(K)
    if fstr(3) == ')';
       x{j+1} = K*(1-0)
    else x{j+1} = integral(K,0,1)
    end
 E(j+1,1) = subs(x{j+1},1);
 j = j+1
 end
end

But the following error occurs,

Undefined operator '*' for input arguments of type 'function_handle'.
Error in solcone1 (line 9)
x{j+1} = K*(1-0);
like image 969
Airapet Avatar asked Jun 08 '18 08:06

Airapet


1 Answers

I am going to ignore the terrible choice of using eval, especially when you can do

x{j+1} = integral(K,0,1);

Read more on why dynamic variables and eval are terrible


Your problem is that matlabFunction is a smartass. When it detects that your function does not have any dependency to x, it gives you a function with empty input arguments @()4.0. As you see, integral does not like that.

A way of solving the problem is detecting this before calling integral. You could check if it has input arguments, and if it does not, then evaluate the integral "by hand"

 ...
 while j < n+1 ;

    K = matlabFunction(subs(po(j,k))) ; 
    fstr=func2str(K);
    if fstr(3)==')'
         x{j+1}=K()*(1-0); % evaluate the integral yourself
    else
         x{j+1} = integral(K,0,1);
    end
    E(j+1,1) = subs(x{j+1});
    j = j+1;
 end
...

The problem is considerably more difficult that I though I was. Either rewrite the entire thing, or using eval:

 ...
 while j < n+1 ;

    K = matlabFunction(subs(po(j,k))) ; 
    fstr=func2str(K);
    if j==0
          eval(sprintf('x%d = K()*(1-0);;',j+1)) ;
    else
          eval(sprintf('x%d = integral(K,0,1);',j+1)) ;
    end
     E(j+1,1) = subs(sprintf('x%d',j+1));
    j = j+1;
 end
...
like image 159
Ander Biguri Avatar answered Sep 19 '22 13:09

Ander Biguri