Theoretical details regarding symbolic integration in Matlab

Dear all,

I am using symbolic integration in Matlab (in particular the "int"-function). I would like to use it as a reference in my research. Therefore, it would be nice to obtain some details which methods Matlab is actually using. Does anybody know? Is there a source to get more details?

Best,

Michael

Answer 1

[Walter Roberson]

The symbolic integration methods are undocumented.

The integration methods are very likely table driven: some pattern matching is done and various integration methods are attempted until eventually one works (or none of them work.)

Forming the grobnar basis is one of the potential steps.

For what it is worth, you can examine code:

regexprep(char(evalin(symengine, 'expose(int)')), '\\n', '\n')

ans = 
    'proc(f, x = null())
       name int;
       local integral, xx, a, b, prop, unevaluated_int, options, optseq, u, errtext;
       option hold;
     begin
       if args(0) = 0 then
         error(message(\"symbolic:int:NoArgument\"))
       end_if; 
       x := context(x); 
       optseq := context(args(3..args(0))); 
       if testtype(x, \"_exprseq\") then
         optseq := (x[2..-1], optseq); 
         x := x[1]
       end_if; 
       options := intlib::getOptions(optseq); 
       if options[Hold] then
         return(hold(int)(context(f), x, optseq))
       end_if; 
       if type(x) = \"_equal\" && type(rhs(x)) = \"_range\" && type(lhs(x)) ~= \"_mult\" then
         xx := op(x, 1); 
         if domtype(xx) ~= DOM_IDENT then
           error(message(\"symbolic:int:VariableMustBeIdentifier2\"))
         end_if; 
         a := op(op(x, 2), 1); 
         b := op(op(x, 2), 2); 
         if protected(xx) ~= None then
           develinfo(1, \"\".xx.\" is protected, can't do assumptions early on\")
         else 
           if options[hold(Assumptions)] ~= FALSE then
             if is(a <= b) = UNKNOWN && intlib::printWarningsFlag && options[NoWarning] ~= TRUE then
               warning(message(\"symbolic:int:Assuming\", a <= b))
             end_if; 
             prop := Dom::Interval([a, b]); 
             if type(prop) = Dom::Interval || prop = R_ then
               if property::hasprop(xx) then
                 if intlib::printWarningsFlag && options[NoWarning] ~= TRUE && getprop(xx) ~= prop then
                   warning(message(\"symbolic:int:ChangingProperties2\", xx, prop, getprop(xx)))
                 end_if; 
                 unassume(xx)
               end_if; 
               if is(b < a) = TRUE then
                 assumeAlso(xx in Dom::Interval([b, a]))
               else 
                 if is(a <= b) = TRUE then
                   assumeAlso(xx in Dom::Interval([a, b]))
                 else 
                   assumeAlso(xx in R_)
                 end_if
               end_if
             end_if
           end_if
         end_if
       else 
         if type(x) = \"_equal\" && type(lhs(x)) = \"_mult\" then
           xx := op(x, 1); 
           a := op(op(x, 2), 1); 
           b := op(op(x, 2), 2); 
           [u, xx] := split(xx, unit::isUnit)[1..2]; 
           return(context(hold(int)(f, xx = a/u..b/u, options))*u)
         else 
           if type(x) = \"_mult\" then
             [u, x] := split(x, unit::isUnit)[1..2]; 
             return(context(hold(int)(f, x, options))*u)
           end_if
         end_if
       end_if; 
       f := hold(hold)(context(f)); 
       if testargs() && has(context([args()]), FAIL) then
         error(message(\"symbolic:int:InvalidArgument\"))
       end_if; 
       if has(context([args()]), undefined) then
         return(undefined)
       end_if; 
       assert(op(f, 0) = hold(hold)); 
       f := op(f); 
       if f::dom::int ~= FAIL then
         return(f::dom::int(f, x, context(args(3..args(0)))))
       else 
         if domtype(f) = DOM_EXPR && domtype(eval(op(f, 0))) = DOM_FUNC_ENV && (eval(op(f, 0)))::int ~= FAIL then
           return((eval(op(f, 0)))::int(f, x, context(args(3..args(0)))))
         else 
           if testtype(f, Type::Arithmetical) = FALSE then
             error(message(\"symbolic:int:InvalidIntegrand\"))
           end_if
         end_if
       end_if; 
       if args(0) = 1 then
         error(message(\"symbolic:int:MissingSecondArgument\"))
       end_if; 
       if testtype(x, \"_exprseq\") then
         options := intlib::getOptions(x[2..-1], options); 
         x := op(x, 1)
       end_if; 
       if domtype(x) = DOM_IDENT then
         integral := hold(int)(f, x, context(args(3..args(0)))); 
         traperror((integral := intlib::int(f, x, options))); 
         integral := subs(integral, hold(intlib::frozenInt) = hold(int), Unsimplified); 
         unevaluated_int := FALSE; 
         if domtype(integral) ~= piecewise || domtype(piecewise::disregardPoints(integral)) ~= piecewise then
           misc::maprec(piecewise::disregardPoints(integral), {\"int\"} = (i -> (if op(i, [2, 1]) = op(DOM_VAR(1,13), 1) then
               DOM_VAR(1,7) := TRUE; 
               misc::breakmap()
             end_if; 
             i)), Unsimplified)
         end_if; 
         if unevaluated_int && ~options[hold(AlwaysExtractIntegrableParts)] = TRUE && 100 + intlib::Simplify::defaultValuation(hold(int)(f, x)) < intlib::Simplify::defaultValuation(integral) then
           return(hold(int)(f, x, intlib::printOptions(options, x)))
         end_if; 
         if ~unevaluated_int && hastype(integral, piecewise) then
           integral := subs(integral, hold(int) = hold(intlib::frozenInt), Unsimplified); 
           integral := misc::maprec(integral, {piecewise} = (p -> piecewise::mapConditions(p, simplify::simplifyCondition))) assumingAlso f in C_; 
           integral := subs(integral, hold(intlib::frozenInt) = hold(int), Unsimplified)
         end_if; 
         return(integral)
       else 
         if type(x) = \"_equal\" && type(op(x, 2)) = \"_range\" then
           intlib::xrngTest(x); 
           integral := hold(int)(f, x, context(args(3..args(0)))); 
           if traperror((integral := intlib::defInt(f, x, options))) = 1028 then
             errtext := getlasterror()[2]; 
             if stringlib::contains(errtext, (message(\"symbolic:int:InvalidBorder\"))::getString()) then
               error(message(\"symbolic:int:InvalidBorder\"))
             end_if; 
             if stringlib::contains(errtext, (message(\"symbolic:int:DifferentRealParts\"))::getString()) then
               error(message(\"symbolic:int:DifferentRealParts\"))
             end_if
           end_if; 
           if ~hastype(integral, \"int\") && ~hastype(integral, \"limit\") then
             integral := intlib::Simplify(integral)
           end_if; 
           integral := subs(integral, hold(intlib::frozenInt) = hold(int), Unsimplified); 
           unevaluated_int := FALSE; 
           if domtype(integral) ~= piecewise then
             misc::maprec(integral, {\"int\"} = (i -> (if op(i, [2, 1]) = op(DOM_VAR(1,13), 1) then
                 DOM_VAR(1,7) := TRUE; 
                 misc::breakmap()
               end_if; 
               i)), Unsimplified)
           end_if; 
           if unevaluated_int && ~options[hold(AlwaysExtractIntegrableParts)] = TRUE && 3*intlib::Simplify::defaultValuation([f, x]) + 200 < intlib::Simplify::defaultValuation(integral) then
             return(hold(int)(f, x, intlib::printOptions(options, x)))
           end_if; 
           return(integral)
         else 
           error(message(\"symbolic:int:InvalidSecondArgument\"))
         end_if
       end_if
     end_proc'