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) neuper@37906: (*new'_c :: "real => real" ("new'_c _" 66)*) neuper@37906: 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: neuper@37906: (*the CAS-command, eg. "Integrate (2*x^^^3, x)"*) neuper@37906: Integrate :: "[real * real] => real" neuper@37906: neuper@37906: (*Script-names*) neuper@37906: IntegrationScript :: "[real,real, real] => real" neuper@37906: ("((Script IntegrationScript (_ _ =))// (_))" 9) neuper@37906: NamedIntegrationScript :: "[real,real, real=>real, bool] => bool" neuper@37906: ("((Script NamedIntegrationScript (_ _ _=))// (_))" 9) 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: neuper@52148: integral_const: "Not (bdv occurs_in u) ==> Integral u D bdv = u * bdv" and neuper@52148: 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: *) neuper@37983: integral_pow: "Integral bdv ^^^ n D bdv = bdv ^^^ (n+1) / (n + 1)" neuper@37906: neuper@37954: ML {* neuper@37972: val thy = @{theory}; neuper@37972: neuper@37954: (** eval functions **) neuper@37954: neuper@37954: val c = Free ("c", HOLogic.realT); neuper@37954: (*.create a new unique variable 'c..' in a term; for use by Calc 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 neuper@37954: | "c"::"_"::is => (case (int_of_str 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 neuper@37954: "c"::"_"::is => (the o int_of_str o implode) is neuper@37954: | _ => 0; neuper@37954: val cs = filter selc (vars term); neuper@37954: in neuper@37954: case cs of neuper@37954: [] => c neuper@37954: | [c] => 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 neuper@37954: (*("new_c", ("Integrate.new'_c", eval_new_c "#new_c_"))*) neuper@37954: fun eval_new_c _ _ (p as (Const ("Integrate.new'_c",_) $ t)) _ = neuper@37954: SOME ((term2str p) ^ " = " ^ term2str (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:*) neuper@37954: (*("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*) neuper@37954: fun eval_add_new_c (_:string) "Integrate.add'_new'_c" p (_:theory) = neuper@37954: let val p' = case p of neuper@41922: Const ("HOL.eq", T) $ lh $ rh => neuper@41922: Const ("HOL.eq", T) $ lh $ mk_add rh (new_c rh) neuper@37954: | p => mk_add p (new_c p) neuper@37954: in SOME ((term2str p) ^ " = " ^ term2str p', neuper@37954: Trueprop $ (mk_equality (p, p'))) neuper@37954: end neuper@37954: | eval_add_new_c _ _ _ _ = NONE; neuper@37954: neuper@37954: neuper@37954: (*("is_f_x", ("Integrate.is'_f'_x", eval_is_f_x "is_f_x_"))*) neuper@37954: fun eval_is_f_x _ _(p as (Const ("Integrate.is'_f'_x", _) neuper@37954: $ arg)) _ = neuper@37954: if is_f_x arg neuper@37954: then SOME ((term2str p) ^ " = True", neuper@48760: Trueprop $ (mk_equality (p, @{term True}))) neuper@37954: else SOME ((term2str p) ^ " = False", neuper@48760: Trueprop $ (mk_equality (p, @{term False}))) neuper@37954: | eval_is_f_x _ _ _ _ = NONE; neuper@37996: *} s1210629013@52145: setup {* KEStore_Elems.add_calcs s1210629013@52145: [("add_new_c", ("Integrate.add'_new'_c", eval_add_new_c "add_new_c_")), s1210629013@52145: ("is_f_x", ("Integrate.is'_f'_x", eval_is_f_x "is_f_idextifier_"))] *} neuper@37996: ML {* neuper@37954: (** rulesets **) neuper@37954: neuper@37954: (*.rulesets for integration.*) neuper@37954: val integration_rules = neuper@37954: Rls {id="integration_rules", preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), neuper@37954: erls = Rls {id="conditions_in_integration_rules", neuper@37954: preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), neuper@37954: erls = Erls, neuper@42451: srls = Erls, calc = [], errpatts = [], neuper@37954: rules = [(*for rewriting conditions in Thm's*) neuper@37954: Calc ("Atools.occurs'_in", neuper@37954: eval_occurs_in "#occurs_in_"), neuper@37969: Thm ("not_true",num_str @{thm not_true}), neuper@37969: Thm ("not_false",@{thm not_false}) neuper@37954: ], neuper@37954: scr = EmptyScr}, neuper@42451: srls = Erls, calc = [], errpatts = [], neuper@37954: rules = [ neuper@37969: Thm ("integral_const",num_str @{thm integral_const}), neuper@37969: Thm ("integral_var",num_str @{thm integral_var}), neuper@37969: Thm ("integral_add",num_str @{thm integral_add}), neuper@37969: Thm ("integral_mult",num_str @{thm integral_mult}), neuper@37969: Thm ("integral_pow",num_str @{thm integral_pow}), neuper@38014: Calc ("Groups.plus_class.plus", eval_binop "#add_")(*for n+1*) neuper@37954: ], neuper@37954: scr = EmptyScr}; neuper@37996: *} neuper@37996: ML {* neuper@37954: val add_new_c = neuper@37954: Seq {id="add_new_c", preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), neuper@37954: erls = Rls {id="conditions_in_add_new_c", neuper@37954: preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), neuper@37954: erls = Erls, neuper@42451: srls = Erls, calc = [], errpatts = [], neuper@37954: rules = [Calc ("Tools.matches", eval_matches""), neuper@37954: Calc ("Integrate.is'_f'_x", neuper@37954: eval_is_f_x "is_f_x_"), neuper@37969: Thm ("not_true",num_str @{thm not_true}), neuper@37996: Thm ("not_false",num_str @{thm not_false}) neuper@37954: ], neuper@37954: scr = EmptyScr}, neuper@42451: srls = Erls, calc = [], errpatts = [], neuper@37969: rules = [ (*Thm ("call_for_new_c", num_str @{thm call_for_new_c}),*) neuper@37954: Cal1 ("Integrate.add'_new'_c", eval_add_new_c "new_c_") neuper@37954: ], neuper@37954: scr = EmptyScr}; neuper@37996: *} neuper@37996: 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 = neuper@37954: Rls {id = "norm_Rational_rls_noadd_fractions", preconds = [], neuper@37954: rew_ord = ("dummy_ord",dummy_ord), neuper@42451: erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [], neuper@52105: rules = [(*Rls_ add_fractions_p_rls,!!!*) neuper@37954: Rls_ (*rat_mult_div_pow original corrected WN051028*) neuper@37954: (Rls {id = "rat_mult_div_pow", preconds = [], neuper@37954: rew_ord = ("dummy_ord",dummy_ord), neuper@37954: erls = (*FIXME.WN051028 e_rls,*) neuper@37954: append_rls "e_rls-is_polyexp" e_rls neuper@37954: [Calc ("Poly.is'_polyexp", neuper@37954: eval_is_polyexp "")], neuper@42451: srls = Erls, calc = [], errpatts = [], neuper@37969: rules = [Thm ("rat_mult",num_str @{thm rat_mult}), neuper@37954: (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*) neuper@37969: Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}), neuper@37954: (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*) neuper@37969: Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r}), neuper@37954: (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*) neuper@37954: neuper@37979: Thm ("real_divide_divide1_mg", neuper@37979: num_str @{thm real_divide_divide1_mg}), neuper@37954: (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*) neuper@37979: Thm ("divide_divide_eq_right", neuper@37979: num_str @{thm divide_divide_eq_right}), neuper@37954: (*"?x / (?y / ?z) = ?x * ?z / ?y"*) neuper@37979: Thm ("divide_divide_eq_left", neuper@37979: num_str @{thm divide_divide_eq_left}), neuper@37954: (*"?x / ?y / ?z = ?x / (?y * ?z)"*) neuper@48789: Calc ("Fields.inverse_class.divide" ,eval_cancel "#divide_e"), neuper@37954: neuper@37969: Thm ("rat_power", num_str @{thm rat_power}) neuper@37954: (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*) neuper@37954: ], neuper@48763: scr = Prog ((term_of o the o (parse thy)) "empty_script") neuper@37954: }), neuper@37954: Rls_ make_rat_poly_with_parentheses, neuper@37954: Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*) neuper@37954: Rls_ rat_reduce_1 neuper@37954: ], neuper@48763: scr = Prog ((term_of o the o (parse thy)) "empty_script") neuper@37954: }:rls; neuper@37954: neuper@37954: (*.for simplify_Integral adapted from 'norm_Rational'.*) neuper@37954: val norm_Rational_noadd_fractions = neuper@37954: Seq {id = "norm_Rational_noadd_fractions", preconds = [], neuper@37954: rew_ord = ("dummy_ord",dummy_ord), neuper@42451: erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [], neuper@37980: rules = [Rls_ discard_minus, neuper@37954: Rls_ rat_mult_poly,(* removes double fractions like a/b/c *) neuper@37954: Rls_ make_rat_poly_with_parentheses, (*WN0510 also in(#)below*) neuper@37954: Rls_ cancel_p_rls, (*FIXME.MG:cancel_p does NOT order sometim*) neuper@37954: Rls_ norm_Rational_rls_noadd_fractions,(* the main rls (#) *) neuper@37979: Rls_ discard_parentheses1 (* mult only *) neuper@37954: ], neuper@48763: scr = Prog ((term_of o the o (parse thy)) "empty_script") neuper@37954: }:rls; 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 = neuper@37954: append_rls "separate_bdv2" neuper@37954: collect_bdv neuper@37969: [Thm ("separate_bdv", num_str @{thm separate_bdv}), neuper@37954: (*"?a * ?bdv / ?b = ?a / ?b * ?bdv"*) neuper@37969: Thm ("separate_bdv_n", num_str @{thm separate_bdv_n}), neuper@37969: Thm ("separate_1_bdv", num_str @{thm separate_1_bdv}), neuper@37954: (*"?bdv / ?b = (1 / ?b) * ?bdv"*) neuper@37969: Thm ("separate_1_bdv_n", num_str @{thm separate_1_bdv_n})(*, neuper@37954: (*"?bdv ^^^ ?n / ?b = 1 / ?b * ?bdv ^^^ ?n"*) neuper@37991: *****Thm ("add_divide_distrib", neuper@37991: *****num_str @{thm add_divide_distrib}) neuper@37954: (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*)----------*) neuper@37954: ]; neuper@37954: val simplify_Integral = neuper@37954: Seq {id = "simplify_Integral", preconds = []:term list, neuper@37954: rew_ord = ("dummy_ord", dummy_ord), neuper@37954: erls = Atools_erls, srls = Erls, neuper@42451: calc = [], errpatts = [], neuper@52062: rules = [Thm ("distrib_right",num_str @{thm distrib_right}), neuper@37954: (*"(?z1.0 + ?z2.0) * ?w = ?z1.0 * ?w + ?z2.0 * ?w"*) neuper@37991: Thm ("add_divide_distrib",num_str @{thm add_divide_distrib}), neuper@37954: (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*) neuper@37954: (*^^^^^ *1* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*) neuper@37954: Rls_ norm_Rational_noadd_fractions, neuper@37954: Rls_ order_add_mult_in, neuper@37954: Rls_ discard_parentheses, neuper@37954: (*Rls_ collect_bdv, from make_polynomial_in*) neuper@37954: Rls_ separate_bdv2, neuper@48789: Calc ("Fields.inverse_class.divide" ,eval_cancel "#divide_e") neuper@37954: ], neuper@37954: scr = EmptyScr}:rls; 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'( neuper@37954: * Seq {id = "", preconds = []:term list, neuper@37954: * rew_ord = ("dummy_ord", dummy_ord), neuper@37954: * erls = Atools_erls, srls = Erls, neuper@37954: * calc = [], (*asm_thm = [],*) neuper@37954: * rules = [Rls_ expand_poly, neuper@37954: * Rls_ order_add_mult_in, neuper@37954: * Rls_ simplify_power, neuper@37954: * Rls_ collect_numerals, neuper@37954: * Rls_ reduce_012, neuper@37969: * Thm ("realpow_oneI",num_str @{thm realpow_oneI}), neuper@37954: * Rls_ discard_parentheses, neuper@37954: * Rls_ collect_bdv, neuper@37954: * (*below inserted from 'make_ratpoly_in'*) neuper@37954: * Rls_ (append_rls "separate_bdv" neuper@37954: * collect_bdv neuper@37969: * [Thm ("separate_bdv", num_str @{thm separate_bdv}), neuper@37954: * (*"?a * ?bdv / ?b = ?a / ?b * ?bdv"*) neuper@37969: * Thm ("separate_bdv_n", num_str @{thm separate_bdv_n}), neuper@37969: * Thm ("separate_1_bdv", num_str @{thm separate_1_bdv}), neuper@37954: * (*"?bdv / ?b = (1 / ?b) * ?bdv"*) neuper@37969: * Thm ("separate_1_bdv_n", num_str @{thm separate_1_bdv_n})(*, neuper@37954: * (*"?bdv ^^^ ?n / ?b = 1 / ?b * ?bdv ^^^ ?n"*) neuper@37991: * Thm ("add_divide_distrib", neuper@37991: * num_str @{thm add_divide_distrib}) neuper@37954: * (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*)*) neuper@37954: * ]), neuper@48789: * Calc ("Fields.inverse_class.divide" ,eval_cancel "#divide_e") neuper@37954: * ], neuper@37954: * scr = EmptyScr neuper@37954: * }:rls); neuper@37954: .......................................................................*) neuper@37954: neuper@37954: val integration = neuper@37954: Seq {id="integration", preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), neuper@37954: erls = Rls {id="conditions_in_integration", neuper@37954: preconds = [], neuper@37954: rew_ord = ("termlessI",termlessI), neuper@37954: erls = Erls, neuper@42451: srls = Erls, calc = [], errpatts = [], neuper@37954: rules = [], neuper@37954: scr = EmptyScr}, neuper@42451: srls = Erls, calc = [], errpatts = [], neuper@37954: rules = [ Rls_ integration_rules, neuper@37954: Rls_ add_new_c, neuper@37954: Rls_ simplify_Integral neuper@37954: ], neuper@37954: scr = EmptyScr}; s1210629013@55444: s1210629013@55444: val prep_rls' = prep_rls @{theory}; neuper@37996: *} neuper@52125: setup {* KEStore_Elems.add_rlss s1210629013@55444: [("integration_rules", (Context.theory_name @{theory}, prep_rls' integration_rules)), s1210629013@55444: ("add_new_c", (Context.theory_name @{theory}, prep_rls' add_new_c)), s1210629013@55444: ("simplify_Integral", (Context.theory_name @{theory}, prep_rls' simplify_Integral)), s1210629013@55444: ("integration", (Context.theory_name @{theory}, prep_rls' integration)), s1210629013@55444: ("separate_bdv2", (Context.theory_name @{theory}, prep_rls' separate_bdv2)), neuper@52125: neuper@52125: ("norm_Rational_noadd_fractions", (Context.theory_name @{theory}, s1210629013@55444: prep_rls' norm_Rational_noadd_fractions)), neuper@52125: ("norm_Rational_rls_noadd_fractions", (Context.theory_name @{theory}, s1210629013@55444: prep_rls' norm_Rational_rls_noadd_fractions))] *} neuper@37954: neuper@37954: (** problems **) s1210629013@55359: setup {* KEStore_Elems.add_pbts s1210629013@55339: [(prep_pbt thy "pbl_fun_integ" [] e_pblID s1210629013@55339: (["integrate","function"], s1210629013@55339: [("#Given" ,["functionTerm f_f", "integrateBy v_v"]), s1210629013@55339: ("#Find" ,["antiDerivative F_F"])], s1210629013@55339: append_rls "e_rls" e_rls [(*for preds in where_*)], s1210629013@55339: SOME "Integrate (f_f, v_v)", s1210629013@55339: [["diff","integration"]])), s1210629013@55339: (*here "named" is used differently from Differentiation"*) s1210629013@55339: (prep_pbt thy "pbl_fun_integ_nam" [] e_pblID s1210629013@55339: (["named","integrate","function"], s1210629013@55339: [("#Given" ,["functionTerm f_f", "integrateBy v_v"]), s1210629013@55339: ("#Find" ,["antiDerivativeName F_F"])], s1210629013@55339: append_rls "e_rls" e_rls [(*for preds in where_*)], s1210629013@55339: SOME "Integrate (f_f, v_v)", s1210629013@55339: [["diff","integration","named"]]))] *} s1210629013@55380: neuper@37954: (** methods **) s1210629013@55373: setup {* KEStore_Elems.add_mets s1210629013@55373: [prep_met thy "met_diffint" [] e_metID s1210629013@55373: (["diff","integration"], s1210629013@55373: [("#Given" ,["functionTerm f_f", "integrateBy v_v"]), ("#Find" ,["antiDerivative F_F"])], s1210629013@55373: {rew_ord'="tless_true", rls'=Atools_erls, calc = [], srls = e_rls, prls=e_rls, s1210629013@55373: crls = Atools_erls, errpats = [], nrls = e_rls}, s1210629013@55373: "Script IntegrationScript (f_f::real) (v_v::real) = " ^ s1210629013@55373: " (let t_t = Take (Integral f_f D v_v) " ^ s1210629013@55373: " in (Rewrite_Set_Inst [(bdv,v_v)] integration False) (t_t::real))"), s1210629013@55373: prep_met thy "met_diffint_named" [] e_metID s1210629013@55373: (["diff","integration","named"], s1210629013@55373: [("#Given" ,["functionTerm f_f", "integrateBy v_v"]), s1210629013@55373: ("#Find" ,["antiDerivativeName F_F"])], s1210629013@55373: {rew_ord'="tless_true", rls'=Atools_erls, calc = [], srls = e_rls, prls=e_rls, s1210629013@55373: crls = Atools_erls, errpats = [], nrls = e_rls}, s1210629013@55373: "Script NamedIntegrationScript (f_f::real) (v_v::real) (F_F::real=>real) = " ^ s1210629013@55373: " (let t_t = Take (F_F v_v = Integral f_f D v_v) " ^ s1210629013@55373: " in ((Try (Rewrite_Set_Inst [(bdv,v_v)] simplify_Integral False)) @@ " ^ s1210629013@55373: " (Rewrite_Set_Inst [(bdv,v_v)] integration False)) t_t) ")] s1210629013@55373: *} neuper@37954: neuper@37906: end