10

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);
Airapet
  • 103
  • 6
  • 3
    It will not fix your problem , but definetly it will improve your programming and future headaches if you woudl just avoid using `eval` completely. Use `x{jj+1}` – Ander Biguri Jun 08 '18 at 10:46

1 Answers1

11

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
...
Ander Biguri
  • 35,140
  • 11
  • 74
  • 120
  • I can't use these brackets in `x{j+1}` unless I specify a cell array. But my `x` isn't defined that way in `po`. – Airapet Jun 08 '18 at 13:27
  • @Airapet A cell array is an array that contains whatever you want on it. In which way is your function specified so it is not within the "whatever you want" definition? `syms x; f=x.^2; g=@(x)(x^2); a{1}=f; a{2}=g;` works fine – Ander Biguri Jun 08 '18 at 13:36
  • When I call the function now I get the error which I posted under EDIT 2 – Airapet Jun 08 '18 at 14:04
  • @Airapet That is because you are a bit bad a copying :P. Its `K()` not `K` – Ander Biguri Jun 08 '18 at 14:20
  • Exactly what I said in the comment.......`x{j+1}=K()*(1-0);` NOT `x{j+1}=K*(1-0);`. `K` is a fucntion, you need to evaluate it, even if it has no inputs. – Ander Biguri Jun 08 '18 at 14:22
  • Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/172771/discussion-between-airapet-and-ander-biguri). – Airapet Jun 08 '18 at 14:28
  • "I am going to ignore the terrible choice of using `eval`... but do read this link because I can't really ignore it at all." LOL! – Cris Luengo Jun 08 '18 at 16:24
  • @Cris it was that ir rewriting the entire thing, which I was la to do :p – Ander Biguri Jun 09 '18 at 10:24