neuper@37906: (* integration over the reals neuper@37906: author: Walther Neuper neuper@37906: 050814, 08:51 neuper@37906: (c) due to copyright terms neuper@37906: *) neuper@37906: neuper@37954: theory Integrate imports Diff begin neuper@37906: neuper@37906: consts neuper@37906: neuper@37906: Integral :: "[real, real]=> real" ("Integral _ D _" 91) walther@60278: (*new_c :: "real => real" ("new_c _" 66)*) walther@60278: is_f_x :: "real => bool" ("_ is'_f'_x" 10) neuper@37906: neuper@37906: (*descriptions in the related problems*) neuper@37996: integrateBy :: "real => una" neuper@37996: antiDerivative :: "real => una" neuper@37996: antiDerivativeName :: "(real => real) => una" neuper@37906: walther@60260: (*the CAS-command, eg. "Integrate (2*x \ 3, x)"*) neuper@37906: Integrate :: "[real * real] => real" neuper@37906: neuper@52148: axiomatization where neuper@37906: (*stated as axioms, todo: prove as theorems neuper@37906: 'bdv' is a constant handled on the meta-level neuper@37906: specifically as a 'bound variable' *) neuper@37906: walther@60269: (*Ambiguous input\<^here> produces 3 parse trees -----------------------------\\*) neuper@52148: integral_const: "Not (bdv occurs_in u) ==> Integral u D bdv = u * bdv" and walther@60242: integral_var: "Integral bdv D bdv = bdv \ 2 / 2" and neuper@37906: neuper@37983: integral_add: "Integral (u + v) D bdv = neuper@52148: (Integral u D bdv) + (Integral v D bdv)" and neuper@37983: integral_mult: "[| Not (bdv occurs_in u); bdv occurs_in v |] ==> neuper@52148: Integral (u * v) D bdv = u * (Integral v D bdv)" and neuper@37906: (*WN080222: this goes into sub-terms, too ... neuper@37983: call_for_new_c: "[| Not (matches (u + new_c v) a); Not (a is_f_x) |] ==> neuper@37954: a = a + new_c a" neuper@37906: *) walther@60242: integral_pow: "Integral bdv \ n D bdv = bdv \ (n+1) / (n + 1)" walther@60269: (*Ambiguous input\<^here> produces 3 parse trees -----------------------------//*) neuper@37906: wneuper@59472: ML \ neuper@37954: (** eval functions **) neuper@37954: neuper@37954: val c = Free ("c", HOLogic.realT); walther@59878: (*.create a new unique variable 'c..' in a term; for use by Rule.Eval in a rls; neuper@37954: an alternative to do this would be '(Try (Calculate new_c_) (new_c es__))' neuper@37954: in the script; this will be possible if currying doesnt take the value neuper@37954: from a variable, but the value '(new_c es__)' itself.*) neuper@37954: fun new_c term = neuper@37954: let fun selc var = neuper@40836: case (Symbol.explode o id_of) var of neuper@37954: "c"::[] => true walther@59875: | "c"::"_"::is => (case (TermC.int_opt_of_string o implode) is of neuper@37954: SOME _ => true neuper@37954: | NONE => false) neuper@37954: | _ => false; neuper@40836: fun get_coeff c = case (Symbol.explode o id_of) c of walther@59875: "c"::"_"::is => (the o TermC.int_opt_of_string o implode) is neuper@37954: | _ => 0; wneuper@59389: val cs = filter selc (TermC.vars term); neuper@37954: in neuper@37954: case cs of neuper@37954: [] => c walther@60269: | [_] => Free ("c_2", HOLogic.realT) neuper@37954: | cs => neuper@37954: let val max_coeff = maxl (map get_coeff cs) neuper@37954: in Free ("c_"^string_of_int (max_coeff + 1), HOLogic.realT) end neuper@37954: end; neuper@37954: neuper@37954: (*WN080222 walther@60278: (*("new_c", ("Integrate.new_c", eval_new_c "#new_c_"))*) walther@60278: fun eval_new_c _ _ (p as (Const ("Integrate.new_c",_) $ t)) _ = walther@59868: SOME ((UnparseC.term p) ^ " = " ^ UnparseC.term (new_c p), neuper@37954: Trueprop $ (mk_equality (p, new_c p))) neuper@37954: | eval_new_c _ _ _ _ = NONE; neuper@37954: *) neuper@37954: neuper@37954: (*WN080222:*) walther@60278: (*("add_new_c", ("Integrate.add_new_c", eval_add_new_c "#add_new_c_")) neuper@37954: add a new c to a term or a fun-equation; neuper@37954: this is _not in_ the term, because only applied to _whole_ term*) walther@60278: fun eval_add_new_c (_:string) "Integrate.add_new_c" p (_:theory) = neuper@37954: let val p' = case p of wenzelm@60309: Const (\<^const_name>\HOL.eq\, T) $ lh $ rh => wenzelm@60309: Const (\<^const_name>\HOL.eq\, T) $ lh $ TermC.mk_add rh (new_c rh) wneuper@59389: | p => TermC.mk_add p (new_c p) walther@59868: in SOME ((UnparseC.term p) ^ " = " ^ UnparseC.term p', wneuper@59390: HOLogic.Trueprop $ (TermC.mk_equality (p, p'))) neuper@37954: end neuper@37954: | eval_add_new_c _ _ _ _ = NONE; neuper@37954: neuper@37954: walther@60278: (*("is_f_x", ("Integrate.is_f_x", eval_is_f_x "is_f_x_"))*) walther@60278: fun eval_is_f_x _ _(p as (Const ("Integrate.is_f_x", _) neuper@37954: $ arg)) _ = wneuper@59389: if TermC.is_f_x arg walther@59868: then SOME ((UnparseC.term p) ^ " = True", wneuper@59390: HOLogic.Trueprop $ (TermC.mk_equality (p, @{term True}))) walther@59868: else SOME ((UnparseC.term p) ^ " = False", wneuper@59390: HOLogic.Trueprop $ (TermC.mk_equality (p, @{term False}))) neuper@37954: | eval_is_f_x _ _ _ _ = NONE; wneuper@59472: \ wneuper@59472: setup \KEStore_Elems.add_calcs wenzelm@60312: [("add_new_c", ("Integrate.add_new_c" (* FIXME proper const!? *), eval_add_new_c "add_new_c_")), wenzelm@60312: ("is_f_x", (\<^const_name>\is_f_x\, eval_is_f_x "is_f_idextifier_"))]\ wneuper@59472: ML \ neuper@37954: (** rulesets **) neuper@37954: neuper@37954: (*.rulesets for integration.*) neuper@37954: val integration_rules = walther@59851: Rule_Def.Repeat {id="integration_rules", preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), walther@59851: erls = Rule_Def.Repeat {id="conditions_in_integration_rules", neuper@37954: preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), walther@59851: erls = Rule_Set.Empty, walther@59851: srls = Rule_Set.Empty, calc = [], errpatts = [], neuper@37954: rules = [(*for rewriting conditions in Thm's*) wenzelm@60294: \<^rule_eval>\Prog_Expr.occurs_in\ (Prog_Expr.eval_occurs_in "#occurs_in_"), wenzelm@60297: \<^rule_thm>\not_true\, wenzelm@60298: \<^rule_thm>\not_false\ neuper@37954: ], walther@59878: scr = Rule.Empty_Prog}, walther@59851: srls = Rule_Set.Empty, calc = [], errpatts = [], neuper@37954: rules = [ wenzelm@60297: \<^rule_thm>\integral_const\, wenzelm@60297: \<^rule_thm>\integral_var\, wenzelm@60297: \<^rule_thm>\integral_add\, wenzelm@60297: \<^rule_thm>\integral_mult\, wenzelm@60297: \<^rule_thm>\integral_pow\, wenzelm@60294: \<^rule_eval>\plus\ (**)(eval_binop "#add_")(*for n+1*) neuper@37954: ], walther@59878: scr = Rule.Empty_Prog}; wneuper@59472: \ wneuper@59472: ML \ neuper@37954: val add_new_c = walther@59878: Rule_Set.Sequence {id="add_new_c", preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), walther@59851: erls = Rule_Def.Repeat {id="conditions_in_add_new_c", neuper@37954: preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), walther@59851: erls = Rule_Set.Empty, walther@59851: srls = Rule_Set.Empty, calc = [], errpatts = [], wenzelm@60294: rules = [\<^rule_eval>\Prog_Expr.matches\ (Prog_Expr.eval_matches""), wenzelm@60294: \<^rule_eval>\Integrate.is_f_x\ (eval_is_f_x "is_f_x_"), wenzelm@60297: \<^rule_thm>\not_true\, wenzelm@60297: \<^rule_thm>\not_false\ neuper@37954: ], walther@59878: scr = Rule.Empty_Prog}, walther@59851: srls = Rule_Set.Empty, calc = [], errpatts = [], wenzelm@60297: rules = [ (*\<^rule_thm>\call_for_new_c\,*) walther@60278: Rule.Cal1 ("Integrate.add_new_c", eval_add_new_c "new_c_") neuper@37954: ], walther@59878: scr = Rule.Empty_Prog}; wneuper@59472: \ wneuper@59472: ML \ neuper@37954: neuper@37954: (*.rulesets for simplifying Integrals.*) neuper@37954: neuper@37954: (*.for simplify_Integral adapted from 'norm_Rational_rls'.*) neuper@37954: val norm_Rational_rls_noadd_fractions = walther@59851: Rule_Def.Repeat {id = "norm_Rational_rls_noadd_fractions", preconds = [], walther@59857: rew_ord = ("dummy_ord",Rewrite_Ord.dummy_ord), walther@59851: erls = norm_rat_erls, srls = Rule_Set.Empty, calc = [], errpatts = [], wneuper@59416: rules = [(*Rule.Rls_ add_fractions_p_rls,!!!*) wneuper@59416: Rule.Rls_ (*rat_mult_div_pow original corrected WN051028*) walther@59851: (Rule_Def.Repeat {id = "rat_mult_div_pow", preconds = [], walther@59857: rew_ord = ("dummy_ord",Rewrite_Ord.dummy_ord), walther@59852: erls = (*FIXME.WN051028 Rule_Set.empty,*) walther@59852: Rule_Set.append_rules "Rule_Set.empty-is_polyexp" Rule_Set.empty wenzelm@60294: [\<^rule_eval>\is_polyexp\ (eval_is_polyexp "")], walther@59851: srls = Rule_Set.Empty, calc = [], errpatts = [], wenzelm@60297: rules = [\<^rule_thm>\rat_mult\, neuper@37954: (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*) wenzelm@60297: \<^rule_thm>\rat_mult_poly_l\, neuper@37954: (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*) wenzelm@60297: \<^rule_thm>\rat_mult_poly_r\, neuper@37954: (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*) neuper@37954: wenzelm@60298: \<^rule_thm>\real_divide_divide1_mg\, neuper@37954: (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*) wenzelm@60298: \<^rule_thm>\divide_divide_eq_right\, neuper@37954: (*"?x / (?y / ?z) = ?x * ?z / ?y"*) wenzelm@60298: \<^rule_thm>\divide_divide_eq_left\, neuper@37954: (*"?x / ?y / ?z = ?x / (?y * ?z)"*) wenzelm@60294: \<^rule_eval>\divide\ (Prog_Expr.eval_cancel "#divide_e"), neuper@37954: wenzelm@60297: \<^rule_thm>\rat_power\ walther@60260: (*"(?a / ?b) \ ?n = ?a \ ?n / ?b \ ?n"*) neuper@37954: ], walther@59878: scr = Rule.Empty_Prog neuper@37954: }), wneuper@59416: Rule.Rls_ make_rat_poly_with_parentheses, wneuper@59416: Rule.Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*) wneuper@59416: Rule.Rls_ rat_reduce_1 neuper@37954: ], walther@59878: scr = Rule.Empty_Prog wneuper@59406: }; neuper@37954: neuper@37954: (*.for simplify_Integral adapted from 'norm_Rational'.*) neuper@37954: val norm_Rational_noadd_fractions = walther@59878: Rule_Set.Sequence {id = "norm_Rational_noadd_fractions", preconds = [], walther@59857: rew_ord = ("dummy_ord",Rewrite_Ord.dummy_ord), walther@59851: erls = norm_rat_erls, srls = Rule_Set.Empty, calc = [], errpatts = [], wneuper@59416: rules = [Rule.Rls_ discard_minus, wneuper@59416: Rule.Rls_ rat_mult_poly,(* removes double fractions like a/b/c *) wneuper@59416: Rule.Rls_ make_rat_poly_with_parentheses, (*WN0510 also in(#)below*) wneuper@59416: Rule.Rls_ cancel_p_rls, (*FIXME.MG:cancel_p does NOT order sometim*) wneuper@59416: Rule.Rls_ norm_Rational_rls_noadd_fractions,(* the main rls (#) *) wneuper@59416: Rule.Rls_ discard_parentheses1 (* mult only *) neuper@37954: ], walther@59878: scr = Rule.Empty_Prog wneuper@59406: }; neuper@37954: neuper@37954: (*.simplify terms before and after Integration such that neuper@37954: ..a.x^2/2 + b.x^3/3.. is made to ..a/2.x^2 + b/3.x^3.. (and NO neuper@37954: common denominator as done by norm_Rational or make_ratpoly_in. neuper@37954: This is a copy from 'make_ratpoly_in' with respective reduction of rules and neuper@37954: *1* expand the term, ie. distribute * and / over + neuper@37954: .*) neuper@37954: val separate_bdv2 = walther@59852: Rule_Set.append_rules "separate_bdv2" neuper@37954: collect_bdv wenzelm@60297: [\<^rule_thm>\separate_bdv\, neuper@37954: (*"?a * ?bdv / ?b = ?a / ?b * ?bdv"*) wenzelm@60297: \<^rule_thm>\separate_bdv_n\, wenzelm@60297: \<^rule_thm>\separate_1_bdv\, neuper@37954: (*"?bdv / ?b = (1 / ?b) * ?bdv"*) wenzelm@60297: \<^rule_thm>\separate_1_bdv_n\(*, walther@60260: (*"?bdv \ ?n / ?b = 1 / ?b * ?bdv \ ?n"*) wenzelm@60298: *****\<^rule_thm>\add_divide_distrib\ neuper@37954: (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*)----------*) neuper@37954: ]; neuper@37954: val simplify_Integral = walther@59878: Rule_Set.Sequence {id = "simplify_Integral", preconds = []:term list, walther@59857: rew_ord = ("dummy_ord", Rewrite_Ord.dummy_ord), walther@59851: erls = Atools_erls, srls = Rule_Set.Empty, neuper@42451: calc = [], errpatts = [], wenzelm@60297: rules = [\<^rule_thm>\distrib_right\, neuper@37954: (*"(?z1.0 + ?z2.0) * ?w = ?z1.0 * ?w + ?z2.0 * ?w"*) wenzelm@60297: \<^rule_thm>\add_divide_distrib\, neuper@37954: (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*) neuper@37954: (*^^^^^ *1* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*) wneuper@59416: Rule.Rls_ norm_Rational_noadd_fractions, wneuper@59416: Rule.Rls_ order_add_mult_in, wneuper@59416: Rule.Rls_ discard_parentheses, wneuper@59416: (*Rule.Rls_ collect_bdv, from make_polynomial_in*) wneuper@59416: Rule.Rls_ separate_bdv2, wenzelm@60294: \<^rule_eval>\divide\ (Prog_Expr.eval_cancel "#divide_e") neuper@37954: ], walther@59878: scr = Rule.Empty_Prog}; neuper@37954: neuper@37954: neuper@37954: (*simplify terms before and after Integration such that neuper@37954: ..a.x^2/2 + b.x^3/3.. is made to ..a/2.x^2 + b/3.x^3.. (and NO neuper@37954: common denominator as done by norm_Rational or make_ratpoly_in. neuper@37954: This is a copy from 'make_polynomial_in' with insertions from neuper@37954: 'make_ratpoly_in' neuper@37954: THIS IS KEPT FOR COMPARISON ............................................ s1210629013@55444: * val simplify_Integral = prep_rls'( walther@59878: * Rule_Set.Sequence {id = "", preconds = []:term list, walther@59857: * rew_ord = ("dummy_ord", Rewrite_Ord.dummy_ord), walther@59851: * erls = Atools_erls, srls = Rule_Set.Empty, neuper@37954: * calc = [], (*asm_thm = [],*) wneuper@59416: * rules = [Rule.Rls_ expand_poly, wneuper@59416: * Rule.Rls_ order_add_mult_in, wneuper@59416: * Rule.Rls_ simplify_power, wneuper@59416: * Rule.Rls_ collect_numerals, wneuper@59416: * Rule.Rls_ reduce_012, wenzelm@60297: * \<^rule_thm>\realpow_oneI\, wneuper@59416: * Rule.Rls_ discard_parentheses, wneuper@59416: * Rule.Rls_ collect_bdv, neuper@37954: * (*below inserted from 'make_ratpoly_in'*) walther@59852: * Rule.Rls_ (Rule_Set.append_rules "separate_bdv" neuper@37954: * collect_bdv wenzelm@60297: * [\<^rule_thm>\separate_bdv\, neuper@37954: * (*"?a * ?bdv / ?b = ?a / ?b * ?bdv"*) wenzelm@60297: * \<^rule_thm>\separate_bdv_n\, wenzelm@60297: * \<^rule_thm>\separate_1_bdv\, neuper@37954: * (*"?bdv / ?b = (1 / ?b) * ?bdv"*) wenzelm@60297: * \<^rule_thm>\separate_1_bdv_n\(*, walther@60260: * (*"?bdv \ ?n / ?b = 1 / ?b * ?bdv \ ?n"*) wenzelm@60298: * \<^rule_thm>\add_divide_distrib\ neuper@37954: * (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*)*) neuper@37954: * ]), wenzelm@60294: * \<^rule_eval>\divide\ (eval_cancel "#divide_e") neuper@37954: * ], walther@59878: * scr = Rule.Empty_Prog wneuper@59406: * }); neuper@37954: .......................................................................*) neuper@37954: neuper@37954: val integration = walther@59878: Rule_Set.Sequence {id="integration", preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), walther@59851: erls = Rule_Def.Repeat {id="conditions_in_integration", neuper@37954: preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), walther@59851: erls = Rule_Set.Empty, walther@59851: srls = Rule_Set.Empty, calc = [], errpatts = [], neuper@37954: rules = [], walther@59878: scr = Rule.Empty_Prog}, walther@59851: srls = Rule_Set.Empty, calc = [], errpatts = [], wneuper@59416: rules = [ Rule.Rls_ integration_rules, wneuper@59416: Rule.Rls_ add_new_c, wneuper@59416: Rule.Rls_ simplify_Integral neuper@37954: ], walther@59878: scr = Rule.Empty_Prog}; s1210629013@55444: walther@59618: val prep_rls' = Auto_Prog.prep_rls @{theory}; wneuper@59472: \ wenzelm@60289: rule_set_knowledge wenzelm@60286: integration_rules = \prep_rls' integration_rules\ and wenzelm@60286: add_new_c = \prep_rls' add_new_c\ and wenzelm@60286: simplify_Integral = \prep_rls' simplify_Integral\ and wenzelm@60286: integration = \prep_rls' integration\ and wenzelm@60286: separate_bdv2 = \prep_rls' separate_bdv2\ and wenzelm@60286: norm_Rational_noadd_fractions = \prep_rls' norm_Rational_noadd_fractions\ and wenzelm@60286: norm_Rational_rls_noadd_fractions = \prep_rls' norm_Rational_rls_noadd_fractions\ neuper@37954: neuper@37954: (** problems **) wenzelm@60306: wenzelm@60306: problem pbl_fun_integ : "integrate/function" = wenzelm@60306: \Rule_Set.append_rules "empty" Rule_Set.empty [(*for preds in where_*)]\ wenzelm@60306: Method: "diff/integration" wenzelm@60306: CAS: "Integrate (f_f, v_v)" wenzelm@60306: Given: "functionTerm f_f" "integrateBy v_v" wenzelm@60306: Find: "antiDerivative F_F" wenzelm@60306: wenzelm@60306: problem pbl_fun_integ_nam : "named/integrate/function" = wenzelm@60306: (*here "named" is used differently from Differentiation"*) wenzelm@60306: \Rule_Set.append_rules "empty" Rule_Set.empty [(*for preds in where_*)]\ wenzelm@60306: Method: "diff/integration/named" wenzelm@60306: CAS: "Integrate (f_f, v_v)" wenzelm@60306: Given: "functionTerm f_f" "integrateBy v_v" wenzelm@60306: Find: "antiDerivativeName F_F" s1210629013@55380: neuper@37954: (** methods **) wneuper@59545: wneuper@59504: partial_function (tailrec) integrate :: "real \ real \ real" wneuper@59504: where walther@59635: "integrate f_f v_v = ( walther@59635: let walther@59635: t_t = Take (Integral f_f D v_v) walther@59635: in walther@59635: (Rewrite_Set_Inst [(''bdv'', v_v)] ''integration'') t_t)" wenzelm@60303: wenzelm@60303: method met_diffint : "diff/integration" = wenzelm@60303: \{rew_ord'="tless_true", rls'=Atools_erls, calc = [], srls = Rule_Set.empty, prls=Rule_Set.empty, wenzelm@60303: crls = Atools_erls, errpats = [], nrls = Rule_Set.empty}\ wenzelm@60303: Program: integrate.simps wenzelm@60303: Given: "functionTerm f_f" "integrateBy v_v" wenzelm@60303: Find: "antiDerivative F_F" wneuper@59545: wneuper@59504: partial_function (tailrec) intergrate_named :: "real \ real \ (real \ real) \ bool" walther@59635: where walther@59635: "intergrate_named f_f v_v F_F =( walther@59635: let walther@59635: t_t = Take (F_F v_v = Integral f_f D v_v) walther@59635: in ( walther@59637: (Try (Rewrite_Set_Inst [(''bdv'', v_v)] ''simplify_Integral'')) #> walther@59635: (Rewrite_Set_Inst [(''bdv'', v_v)] ''integration'') walther@59635: ) t_t)" wenzelm@60303: wenzelm@60303: method met_diffint_named : "diff/integration/named" = wenzelm@60303: \{rew_ord'="tless_true", rls'=Atools_erls, calc = [], srls = Rule_Set.empty, prls=Rule_Set.empty, wenzelm@60303: crls = Atools_erls, errpats = [], nrls = Rule_Set.empty}\ wenzelm@60303: Program: intergrate_named.simps wenzelm@60303: Given: "functionTerm f_f" "integrateBy v_v" wenzelm@60303: Find: "antiDerivativeName F_F" wenzelm@60303: wenzelm@60303: ML \ walther@60278: \ ML \ wneuper@59472: \ neuper@37954: neuper@37906: end