haftmann@28952: (* Author : Jacques D. Fleuriot paulson@12224: Copyright : 2001 University of Edinburgh paulson@15079: Conversion to Isar and new proofs by Lawrence C Paulson, 2004 bulwahn@41368: Conversion of Mac Laurin to Isar by Lukas Bulwahn and Bernhard Häupler, 2005 paulson@12224: *) paulson@12224: paulson@15944: header{*MacLaurin Series*} paulson@15944: nipkow@15131: theory MacLaurin chaieb@29748: imports Transcendental nipkow@15131: begin obua@14738: paulson@15079: subsection{*Maclaurin's Theorem with Lagrange Form of Remainder*} paulson@15079: paulson@15079: text{*This is a very long, messy proof even now that it's been broken down paulson@15079: into lemmas.*} paulson@15079: paulson@15079: lemma Maclaurin_lemma: paulson@15079: "0 < h ==> nipkow@15539: \B. f h = (\m=0.. fact (Suc m - n) = (Suc m - n) * fact (m - n)" avigad@32031: by (subst fact_reduce_nat, auto) avigad@32031: paulson@15079: lemma Maclaurin_lemma2: hoelzl@41412: fixes B bulwahn@41368: assumes DERIV : "\m t. m < n \ 0\t \ t\h \ DERIV (diff m) t :> diff (Suc m) t" hoelzl@41412: and INIT : "n = Suc k" hoelzl@41412: defines "difg \ (\m t. diff m t - ((\p = 0.. (\m t. diff m t - ?difg m t)") bulwahn@41368: shows "\m t. m < n & 0 \ t & t \ h --> DERIV (difg m) t :> difg (Suc m) t" hoelzl@41412: proof (rule allI impI)+ hoelzl@41412: fix m t assume INIT2: "m < n & 0 \ t & t \ h" hoelzl@41412: have "DERIV (difg m) t :> diff (Suc m) t - hoelzl@41412: ((\x = 0..x = 0..x = 0..x::nat. real (fact x) + real x * real (fact x) \ 0" hoelzl@41412: by (metis fact_gt_zero_nat not_add_less1 real_of_nat_add real_of_nat_mult real_of_nat_zero_iff) hoelzl@41412: have "\x. real (Suc x) * t ^ x * diff (Suc m + x) 0 / real (fact (Suc x)) = hoelzl@41412: diff (Suc m + x) 0 * t^x / real (fact x)" hoelzl@41412: by (auto simp: field_simps real_of_nat_Suc fact_neq_0 intro!: nonzero_divide_eq_eq[THEN iffD2]) hoelzl@41412: moreover hoelzl@41412: have "real (n - m) * t ^ (n - Suc m) * B / real (fact (n - m)) = hoelzl@41412: B * (t ^ (n - Suc m) / real (fact (n - Suc m)))" hoelzl@41412: using `0 < n - m` by (simp add: fact_reduce_nat) hoelzl@41412: ultimately show "DERIV (difg m) t :> difg (Suc m) t" hoelzl@41412: unfolding difg_def by simp bulwahn@41368: qed avigad@32031: paulson@15079: lemma Maclaurin: huffman@29187: assumes h: "0 < h" huffman@29187: assumes n: "0 < n" huffman@29187: assumes diff_0: "diff 0 = f" huffman@29187: assumes diff_Suc: huffman@29187: "\m t. m < n & 0 \ t & t \ h --> DERIV (diff m) t :> diff (Suc m) t" huffman@29187: shows huffman@29187: "\t. 0 < t & t < h & paulson@15079: f h = nipkow@15539: setsum (%m. (diff m 0 / real (fact m)) * h ^ m) {0..m = 0..real) / real (fact m) * h ^ m) + huffman@29187: B * (h ^ n / real (fact n))" huffman@29187: using Maclaurin_lemma [OF h] .. huffman@29187: hoelzl@41412: def g \ "(\t. f t - hoelzl@41412: (setsum (\m. (diff m 0 / real(fact m)) * t^m) {0.. "(%m t. diff m t - huffman@29187: (setsum (%p. (diff (m + p) 0 / real (fact p)) * (t ^ p)) {0..(m\nat) t\real. huffman@29187: m < n \ (0\real) \ t \ t \ h \ DERIV (difg m) t :> difg (Suc m) t" hoelzl@41412: using diff_Suc m unfolding difg_def by (rule Maclaurin_lemma2) huffman@29187: huffman@29187: have difg_eq_0: "\m. m < n --> difg m 0 = 0" huffman@29187: apply clarify huffman@29187: apply (simp add: m difg_def) huffman@29187: apply (frule less_iff_Suc_add [THEN iffD1], clarify) huffman@29187: apply (simp del: setsum_op_ivl_Suc) huffman@30019: apply (insert sumr_offset4 [of "Suc 0"]) avigad@32039: apply (simp del: setsum_op_ivl_Suc fact_Suc) huffman@29187: done huffman@29187: huffman@29187: have isCont_difg: "\m x. \m < n; 0 \ x; x \ h\ \ isCont (difg m) x" huffman@29187: by (rule DERIV_isCont [OF difg_Suc [rule_format]]) simp huffman@29187: huffman@29187: have differentiable_difg: huffman@29187: "\m x. \m < n; 0 \ x; x \ h\ \ difg m differentiable x" huffman@29187: by (rule differentiableI [OF difg_Suc [rule_format]]) simp huffman@29187: huffman@29187: have difg_Suc_eq_0: "\m t. \m < n; 0 \ t; t \ h; DERIV (difg m) t :> 0\ huffman@29187: \ difg (Suc m) t = 0" huffman@29187: by (rule DERIV_unique [OF difg_Suc [rule_format]]) simp huffman@29187: huffman@29187: have "m < n" using m by simp huffman@29187: huffman@29187: have "\t. 0 < t \ t < h \ DERIV (difg m) t :> 0" huffman@29187: using `m < n` huffman@29187: proof (induct m) hoelzl@41412: case 0 huffman@29187: show ?case huffman@29187: proof (rule Rolle) huffman@29187: show "0 < h" by fact huffman@29187: show "difg 0 0 = difg 0 h" by (simp add: difg_0 g2) huffman@29187: show "\x. 0 \ x \ x \ h \ isCont (difg (0\nat)) x" huffman@29187: by (simp add: isCont_difg n) huffman@29187: show "\x. 0 < x \ x < h \ difg (0\nat) differentiable x" huffman@29187: by (simp add: differentiable_difg n) huffman@29187: qed huffman@29187: next hoelzl@41412: case (Suc m') huffman@29187: hence "\t. 0 < t \ t < h \ DERIV (difg m') t :> 0" by simp huffman@29187: then obtain t where t: "0 < t" "t < h" "DERIV (difg m') t :> 0" by fast huffman@29187: have "\t'. 0 < t' \ t' < t \ DERIV (difg (Suc m')) t' :> 0" huffman@29187: proof (rule Rolle) huffman@29187: show "0 < t" by fact huffman@29187: show "difg (Suc m') 0 = difg (Suc m') t" huffman@29187: using t `Suc m' < n` by (simp add: difg_Suc_eq_0 difg_eq_0) huffman@29187: show "\x. 0 \ x \ x \ t \ isCont (difg (Suc m')) x" huffman@29187: using `t < h` `Suc m' < n` by (simp add: isCont_difg) huffman@29187: show "\x. 0 < x \ x < t \ difg (Suc m') differentiable x" huffman@29187: using `t < h` `Suc m' < n` by (simp add: differentiable_difg) huffman@29187: qed huffman@29187: thus ?case huffman@29187: using `t < h` by auto huffman@29187: qed huffman@29187: huffman@29187: then obtain t where "0 < t" "t < h" "DERIV (difg m) t :> 0" by fast huffman@29187: huffman@29187: hence "difg (Suc m) t = 0" huffman@29187: using `m < n` by (simp add: difg_Suc_eq_0) huffman@29187: huffman@29187: show ?thesis huffman@29187: proof (intro exI conjI) huffman@29187: show "0 < t" by fact huffman@29187: show "t < h" by fact huffman@29187: show "f h = huffman@29187: (\m = 0..0 & diff 0 = f & nipkow@25134: (\m t. m < n & 0 \ t & t \ h --> DERIV (diff m) t :> diff (Suc m) t) nipkow@25134: --> (\t. 0 < t & t < h & nipkow@25134: f h = (\m=0..m t. bulwahn@41368: m < n & 0 \ t & t \ h --> DERIV (diff m) t :> diff (Suc m) t" bulwahn@41368: shows "\t. 0 < t \ t \ h \ f h = bulwahn@41368: (\m=0.. 0" by simp bulwahn@41368: from INIT1 this INIT2 DERIV have "\t>0. t < h \ bulwahn@41368: f h = hoelzl@41412: (\m = 0..m t. paulson@15079: m < n & 0 \ t & t \ h --> DERIV (diff m) t :> diff (Suc m) t) paulson@15079: --> (\t. 0 < t & paulson@15079: t \ h & paulson@15079: f h = nipkow@15539: (\m=0..m t. m < n & h \ t & t \ 0 --> DERIV (diff m) t :> diff (Suc m) t" bulwahn@41368: shows "\t. h < t & t < 0 & bulwahn@41368: f h = (\m=0..t>0. t < - h \ bulwahn@41368: f (- (- h)) = bulwahn@41368: (\m = 0.. bulwahn@41368: - t < 0 \ bulwahn@41368: f h = bulwahn@41368: (\m = 0.. 0 & diff 0 = f & paulson@15079: (\m t. paulson@15079: m < n & h \ t & t \ 0 --> DERIV (diff m) t :> diff (Suc m) t)) paulson@15079: --> (\t. h < t & paulson@15079: t < 0 & paulson@15079: f h = nipkow@15539: (\m=0..0 \ nipkow@25134: diff 0 0 = nipkow@25134: (\m = 0..m t. m < n & abs t \ abs x --> DERIV (diff m) t :> diff (Suc m) t" bulwahn@41368: shows "\t. abs t \ abs x & paulson@15079: f x = nipkow@15539: (\m=0..t. _ \ f x = ?f x t") hoelzl@41412: proof cases hoelzl@41412: assume "n = 0" with `diff 0 = f` show ?thesis by force bulwahn@41368: next hoelzl@41412: assume "n \ 0" hoelzl@41412: show ?thesis hoelzl@41412: proof (cases rule: linorder_cases) hoelzl@41412: assume "x = 0" with `n \ 0` `diff 0 = f` DERIV hoelzl@41412: have "\0\ \ \x\ \ f x = ?f x 0" by (force simp add: Maclaurin_bi_le_lemma) hoelzl@41412: thus ?thesis .. bulwahn@41368: next hoelzl@41412: assume "x < 0" hoelzl@41412: with `n \ 0` DERIV hoelzl@41412: have "\t>x. t < 0 \ diff 0 x = ?f x t" by (intro Maclaurin_minus) auto hoelzl@41412: then guess t .. hoelzl@41412: with `x < 0` `diff 0 = f` have "\t\ \ \x\ \ f x = ?f x t" by simp hoelzl@41412: thus ?thesis .. hoelzl@41412: next hoelzl@41412: assume "x > 0" hoelzl@41412: with `n \ 0` `diff 0 = f` DERIV hoelzl@41412: have "\t>0. t < x \ diff 0 x = ?f x t" by (intro Maclaurin) auto hoelzl@41412: then guess t .. hoelzl@41412: with `x > 0` `diff 0 = f` have "\t\ \ \x\ \ f x = ?f x t" by simp hoelzl@41412: thus ?thesis .. bulwahn@41368: qed bulwahn@41368: qed bulwahn@41368: paulson@15079: lemma Maclaurin_all_lt: bulwahn@41368: assumes INIT1: "diff 0 = f" and INIT2: "0 < n" and INIT3: "x \ 0" bulwahn@41368: and DERIV: "\m x. DERIV (diff m) x :> diff(Suc m) x" hoelzl@41412: shows "\t. 0 < abs t & abs t < abs x & f x = hoelzl@41412: (\m=0..t. _ \ _ \ f x = ?f x t") hoelzl@41412: proof (cases rule: linorder_cases) hoelzl@41412: assume "x = 0" with INIT3 show "?thesis".. hoelzl@41412: next hoelzl@41412: assume "x < 0" hoelzl@41412: with assms have "\t>x. t < 0 \ f x = ?f x t" by (intro Maclaurin_minus) auto hoelzl@41412: then guess t .. hoelzl@41412: with `x < 0` have "0 < \t\ \ \t\ < \x\ \ f x = ?f x t" by simp hoelzl@41412: thus ?thesis .. hoelzl@41412: next hoelzl@41412: assume "x > 0" hoelzl@41412: with assms have "\t>0. t < x \ f x = ?f x t " by (intro Maclaurin) auto hoelzl@41412: then guess t .. hoelzl@41412: with `x > 0` have "0 < \t\ \ \t\ < \x\ \ f x = ?f x t" by simp hoelzl@41412: thus ?thesis .. bulwahn@41368: qed bulwahn@41368: paulson@15079: paulson@15079: lemma Maclaurin_all_lt_objl: paulson@15079: "diff 0 = f & paulson@15079: (\m x. DERIV (diff m) x :> diff(Suc m) x) & nipkow@25162: x ~= 0 & n > 0 paulson@15079: --> (\t. 0 < abs t & abs t < abs x & nipkow@15539: f x = (\m=0.. n \ 0 --> nipkow@15539: (\m=0..m x. DERIV (diff m) x :> diff (Suc m) x" hoelzl@41412: shows "\t. abs t \ abs x & f x = hoelzl@41412: (\m=0..t. _ \ f x = ?f x t") hoelzl@41412: proof cases hoelzl@41412: assume "n = 0" with INIT show ?thesis by force bulwahn@41368: next hoelzl@41412: assume "n \ 0" hoelzl@41412: show ?thesis hoelzl@41412: proof cases hoelzl@41412: assume "x = 0" hoelzl@41412: with `n \ 0` have "(\m = 0.. 0` have " \0\ \ \x\ \ f x = ?f x 0" by force hoelzl@41412: thus ?thesis .. hoelzl@41412: next hoelzl@41412: assume "x \ 0" hoelzl@41412: with INIT `n \ 0` DERIV have "\t. 0 < \t\ \ \t\ < \x\ \ f x = ?f x t" hoelzl@41412: by (intro Maclaurin_all_lt) auto hoelzl@41412: then guess t .. hoelzl@41412: hence "\t\ \ \x\ \ f x = ?f x t" by simp hoelzl@41412: thus ?thesis .. bulwahn@41368: qed bulwahn@41368: qed bulwahn@41368: paulson@15079: lemma Maclaurin_all_le_objl: "diff 0 = f & paulson@15079: (\m x. DERIV (diff m) x :> diff (Suc m) x) paulson@15079: --> (\t. abs t \ abs x & nipkow@15539: f x = (\m=0.. 0 |] paulson@15079: ==> (\t. 0 < abs t & paulson@15079: abs t < abs x & nipkow@15539: exp x = (\m=0..t. abs t \ abs x & nipkow@15539: exp x = (\m=0..0 --> Suc (Suc (2 * n - 2)) = 2*n" paulson@15251: by (induct "n", auto) paulson@15079: paulson@15079: lemma lemma_Suc_Suc_4n_diff_2 [rule_format, simp]: nipkow@25134: "n\0 --> Suc (Suc (4*n - 2)) = 4*n" paulson@15251: by (induct "n", auto) paulson@15079: paulson@15079: lemma Suc_mult_two_diff_one [rule_format, simp]: nipkow@25134: "n\0 --> Suc (2 * n - 1) = 2*n" paulson@15251: by (induct "n", auto) paulson@15079: paulson@15234: paulson@15234: text{*It is unclear why so many variant results are needed.*} paulson@15079: huffman@36974: lemma sin_expansion_lemma: hoelzl@41412: "sin (x + real (Suc m) * pi / 2) = huffman@36974: cos (x + real (m) * pi / 2)" huffman@36974: by (simp only: cos_add sin_add real_of_nat_Suc add_divide_distrib left_distrib, auto) huffman@36974: huffman@45166: lemma sin_coeff_0 [simp]: "sin_coeff 0 = 0" huffman@45166: unfolding sin_coeff_def by simp (* TODO: move *) huffman@45166: paulson@15079: lemma Maclaurin_sin_expansion2: paulson@15079: "\t. abs t \ abs x & paulson@15079: sin x = huffman@45166: (\m=0..t. sin x = huffman@45166: (\m=0.. 0; 0 < x |] ==> paulson@15079: \t. 0 < t & t < x & paulson@15079: sin x = huffman@45166: (\m=0.. paulson@15079: \t. 0 < t & t \ x & paulson@15079: sin x = huffman@45166: (\m=0..m=0..<(Suc n). cos_coeff m * 0 ^ m) = 1" paulson@15251: by (induct "n", auto) paulson@15079: huffman@36974: lemma cos_expansion_lemma: huffman@36974: "cos (x + real(Suc m) * pi / 2) = -sin (x + real m * pi / 2)" huffman@36974: by (simp only: cos_add sin_add real_of_nat_Suc left_distrib add_divide_distrib, auto) huffman@36974: paulson@15079: lemma Maclaurin_cos_expansion: paulson@15079: "\t. abs t \ abs x & paulson@15079: cos x = huffman@45166: (\m=0.. 0 |] ==> paulson@15079: \t. 0 < t & t < x & paulson@15079: cos x = huffman@45166: (\m=0.. 0 |] ==> paulson@15079: \t. x < t & t < 0 & paulson@15079: cos x = huffman@45166: (\m=0.. (v::real) |] ==> \(x + u) - y\ \ v" paulson@15079: by auto paulson@15079: paulson@15079: lemma Maclaurin_sin_bound: huffman@45166: "abs(sin x - (\m=0.. inverse(real (fact n)) * \x\ ^ n" obua@14738: proof - paulson@15079: have "!! x (y::real). x \ 1 \ 0 \ y \ x * y \ 1 * y" obua@14738: by (rule_tac mult_right_mono,simp_all) obua@14738: note est = this[simplified] huffman@22985: let ?diff = "\(n::nat) x. if n mod 4 = 0 then sin(x) else if n mod 4 = 1 then cos(x) else if n mod 4 = 2 then -sin(x) else -cos(x)" huffman@22985: have diff_0: "?diff 0 = sin" by simp huffman@22985: have DERIV_diff: "\m x. DERIV (?diff m) x :> ?diff (Suc m) x" huffman@22985: apply (clarify) huffman@22985: apply (subst (1 2 3) mod_Suc_eq_Suc_mod) huffman@22985: apply (cut_tac m=m in mod_exhaust_less_4) hoelzl@31880: apply (safe, auto intro!: DERIV_intros) huffman@22985: done huffman@22985: from Maclaurin_all_le [OF diff_0 DERIV_diff] huffman@22985: obtain t where t1: "\t\ \ \x\" and huffman@22985: t2: "sin x = (\m = 0..m. ?diff m 0 = (if even m then 0 huffman@23177: else -1 ^ ((m - Suc 0) div 2))" huffman@22985: apply (subst even_even_mod_4_iff) huffman@22985: apply (cut_tac m=m in mod_exhaust_less_4) huffman@22985: apply (elim disjE, simp_all) huffman@22985: apply (safe dest!: mod_eqD, simp_all) huffman@22985: done obua@14738: show ?thesis huffman@45166: unfolding sin_coeff_def huffman@22985: apply (subst t2) paulson@15079: apply (rule sin_bound_lemma) nipkow@15536: apply (rule setsum_cong[OF refl]) huffman@22985: apply (subst diff_m_0, simp) paulson@15079: apply (auto intro: mult_right_mono [where b=1, simplified] mult_right_mono hoelzl@41412: simp add: est mult_nonneg_nonneg mult_ac divide_inverse paulson@16924: power_abs [symmetric] abs_mult) obua@14738: done obua@14738: qed obua@14738: paulson@15079: end