doc-src/isac/jrocnik/Test_Z_Transform.thy
author Jan Rocnik <jan.rocnik@student.tugraz.at>
Tue, 06 Sep 2011 15:57:25 +0200
branchdecompose-isar
changeset 42245 3d1f27a0f8a4
parent 42244 e8be2cec4ec5
child 42247 56eebac3d65d
permissions -rwxr-xr-x
tuned z-transform-test, add fourier, tuned bakkarb.
     1 theory Test_Z_Transform imports Isac begin
     2 
     3 section {*trials towards Z transform *}
     4 subsection {*terms*}
     5 ML {*
     6 @{term "1 < || z ||"};
     7 @{term "z / (z - 1)"};
     8 @{term "-u -n - 1"};
     9 @{term "-u [-n - 1]"}; (*[ ] denotes lists !!!*)
    10 @{term "z /(z - 1) = -u [-n - 1]"};
    11 @{term "1 < || z || ==> z / (z - 1) = -u [-n - 1]"};
    12 term2str @{term "1 < || z || ==> z / (z - 1) = -u [-n - 1]"};
    13 *}
    14 ML {*
    15 (*alpha -->  "</alpha>" *)
    16 
    17 @{term "\<alpha> "};
    18 @{term "\<delta> "};
    19 @{term "\<phi> "};
    20 @{term "\<rho> "};
    21 term2str @{term "\<rho> "};
    22 *}
    23 
    24 subsection {*rules*}
    25 (*axiomatization "z / (z - 1) = -u [-n - 1]" Illegal variable name: "z / (z - 1) = -u [-n - 1]" *)
    26 (*definition     "z / (z - 1) = -u [-n - 1]" Bad head of lhs: existing constant "op /"*)
    27 axiomatization where 
    28   rule1: "1 = \<delta>[n]" and
    29   rule2: "|| z || > 1 ==> z / (z - 1) = u [n]" and
    30   rule3: "|| z || < 1 ==> z / (z - 1) = -u [-n - 1]" and 
    31   rule4: "|| z || > || \<alpha> || ==> z / (z - \<alpha>) = \<alpha>^n * u [n]" and
    32   rule5: "|| z || < || \<alpha> || ==> z / (z - \<alpha>) = -(\<alpha>^n) * u [-n - 1]" and
    33   rule6: "|| z || > 1 ==> z/(z - 1)^2 = n * u [n]"
    34 ML {*
    35 @{thm rule1};
    36 @{thm rule2};
    37 @{thm rule3};
    38 @{thm rule4};
    39 *}
    40 
    41 subsection {*apply rules*}
    42 ML {*
    43 val inverse_Z = append_rls "inverse_Z" e_rls
    44   [ Thm  ("rule3",num_str @{thm rule3}),
    45     Thm  ("rule4",num_str @{thm rule4}),
    46     Thm  ("rule1",num_str @{thm rule1})   
    47   ];
    48 
    49 val t = str2term "z / (z - 1) + z / (z - \<alpha>) + 1";
    50 val SOME (t', asm) = rewrite_set_ thy true inverse_Z t;
    51 term2str t' = "z / (z - ?\<delta> [?n]) + z / (z - \<alpha>) + ?\<delta> [?n]"; (*attention rule1 !!!*)
    52 *}
    53 ML {*
    54 val (thy, ro, er) = (@{theory}, tless_true, eval_rls);
    55 *}
    56 ML {*
    57 val SOME (t, asm1) = rewrite_ thy ro er true (num_str @{thm rule3}) t;
    58 term2str t = "- ?u [- ?n - 1] + z / (z - \<alpha>) + 1"; (*- real *)
    59 term2str t;
    60 *}
    61 ML {*
    62 val SOME (t, asm2) = rewrite_ thy ro er true (num_str @{thm rule4}) t;
    63 term2str t = "- ?u [- ?n - 1] + \<alpha> ^ ?n * ?u [?n] + 1"; (*- real *)
    64 term2str t;
    65 *}
    66 ML {*
    67 val SOME (t, asm3) = rewrite_ thy ro er true (num_str @{thm rule1}) t;
    68 term2str t = "- ?u [- ?n - 1] + \<alpha> ^ ?n * ?u [?n] + ?\<delta> [?n]"; (*- real *)
    69 term2str t;
    70 *}
    71 ML {*
    72 terms2str (asm1 @ asm2 @ asm3);
    73 *}
    74 
    75 subsection {*prepare expression*}
    76 ML {*
    77 val ctxt = ProofContext.init_global @{theory Isac};
    78 val ctxt = declare_constraints' [@{term "z::real"}] ctxt;
    79 val SOME expression = parseNEW ctxt "3 / (-1/8 + -1/4*z + z^2)";
    80 term2str expression;
    81 *}
    82 
    83 subsection {*solve equation*}
    84 ML {*(*from test/Tools/isac/Minisubpbl/100-init-rootpbl.sml*)
    85 "----------- Minisubplb/100-init-rootpbl.sml ---------------------";
    86 val denominator = parseNEW ctxt "z^2 - 1/4*z - 1/8 = 0";
    87 val fmz = ["equality (z^2 - 1/4*z - 1/8 = (0::real))", "solveFor z","solutions L"];
    88 val (dI',pI',mI') =("Isac", ["univariate","equation"], ["no_met"]);
    89 *}
    90 ML {*
    91 val (p,_,f,nxt,_,pt) = CalcTreeTEST [(fmz, (dI',pI',mI'))];
    92 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
    93 (*[
    94 (([], Frm), solve (z ^ 2 - 1 / 4 * z - 1 / 8 = 0, z)),
    95 (([1], Frm), z ^ 2 - 1 / 4 * z - 1 / 8 = 0),              bad rewrite order
    96 (([1], Res), -1 / 8 + z ^ 2 + -1 / 4 * z = 0),            bad pattern
    97 (([2], Pbl), solve (-1 / 8 + z ^ 2 + -1 / 4 * z = 0, z)),
    98 (([2,1], Pbl), solve (-1 / 8 + z ^ 2 + -1 / 4 * z = 0, z)),
    99 (([2,1,1], Pbl), solve (-1 / 8 + z ^ 2 + -1 / 4 * z = 0, z)),
   100 (([2,1,1,1], Frm), -1 / 8 + z ^ 2 + -1 / 4 * z = 0)] 
   101 *)
   102 *}
   103 ML {*
   104 val denominator = parseNEW ctxt "-1/8 + -1/4*z + z^2 = 0";
   105 (*ergebnis: [gleichung, was tun?, lösung]*)
   106 val fmz = ["equality (-1/8 + -1/4*z + z^2 = (0::real))", "solveFor z","solutions L"];
   107 (*liste der theoreme die zum lösen benötigt werden, aus isac, keine spezielle methode (no met)*)
   108 val (dI',pI',mI') =
   109   ("Isac", ["pqFormula","degree_2","polynomial","univariate","equation"], ["no_met"]);
   110 (*schritte abarbeiten*)
   111 val (p,_,f,nxt,_,pt) = CalcTreeTEST [(fmz, (dI',pI',mI'))];
   112 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   113 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   114 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   115 val (p,_,f,nxt,_,pt) = me nxt p [] pt; (*val nxt = ("Empty_Tac", ...): tac'_*)
   116 show_pt pt;
   117 *}
   118 
   119 subsection {*partial fraction decomposition*}
   120 subsubsection {*solution of the equation*}
   121 ML {*
   122 val SOME solutions = parseNEW ctxt "[z=1/2, z=-1/4]";
   123 term2str solutions;
   124 atomty solutions;
   125 *}
   126 
   127 subsubsection {*get solutions out of list*}
   128 text {*in isac's CTP-based programming language: let s_1 = NTH 1 solutions; s_2 = NTH 2...*}
   129 ML {*
   130 val Const ("List.list.Cons", _) $ s_1 $ (Const ("List.list.Cons", _) $
   131       s_2 $ Const ("List.list.Nil", _)) = solutions;
   132 term2str s_1;
   133 term2str s_2;
   134 *}
   135 
   136 ML {* (*Solutions as Denominator --> Denominator1 = z - Zeropoint1, Denominator2 = z-Zeropoint2,...*)
   137 val xx = HOLogic.dest_eq s_1;
   138 val s_1' = HOLogic.mk_binop "Groups.minus_class.minus" xx;
   139 val xx = HOLogic.dest_eq s_2;
   140 val s_2' = HOLogic.mk_binop "Groups.minus_class.minus" xx;
   141 term2str s_1';
   142 term2str s_2';
   143 *}
   144 
   145 subsubsection {*build expression*}
   146 text {*in isac's CTP-based programming language: let s_1 = Take numerator / (s_1 * s_2)*}
   147 ML {*
   148 (*The Main Denominator is the multiplikation of the partial fraction denominators*)
   149 val denominator' = HOLogic.mk_binop "Groups.times_class.times" (s_1', s_2') ;
   150 val SOME numerator = parseNEW ctxt "3::real";
   151 
   152 val expression' = HOLogic.mk_binop "Rings.inverse_class.divide" (numerator, denominator');
   153 term2str expression';
   154 *}
   155 
   156 subsubsection {*Ansatz - create partial fractions out of our expression*}
   157 
   158 axiomatization where
   159   ansatz2: "n / (a*b) = A/a + B/(b::real)" and
   160   multiply_eq2: "(n / (a*b) = A/a + B/b) = (a*b*(n  / (a*b)) = a*b*(A/a + B/b))"
   161 
   162 ML {*
   163 (*we use our ansatz2 to rewrite our expression and get an equilation with our expression on the left and the partial fractions of it on the right side*)
   164 val SOME (t1,_) = rewrite_ @{theory Isac} e_rew_ord e_rls false @{thm ansatz2} expression';
   165 term2str t1;
   166 atomty t1;
   167 val eq1 = HOLogic.mk_eq (expression', t1);
   168 term2str eq1;
   169 *}
   170 ML {*
   171 (*eliminate the demoninators by multiplying the left and the right side with the main denominator*)
   172 val SOME (eq2,_) = rewrite_ @{theory Isac} e_rew_ord e_rls false @{thm multiply_eq2} eq1;
   173 term2str eq2;
   174 *}
   175 ML {*
   176 (*simplificatoin*)
   177 val SOME (eq3,_) = rewrite_set_ @{theory Isac} false norm_Rational eq2;
   178 term2str eq3; (*?A ?B not simplified*)
   179 *}
   180 ML {*
   181 val SOME fract1 =
   182   parseNEW ctxt "(z - 1 / 2) * (z - -1 / 4) * (A / (z - 1 / 2) + B / (z - -1 / 4))"; (*A B !*)
   183 val SOME (fract2,_) = rewrite_set_ @{theory Isac} false norm_Rational fract1;
   184 term2str fract2 = "(A + -2 * B + 4 * A * z + 4 * B * z) / 4";
   185 (*term2str fract2 = "A * (1 / 4 + z) + B * (-1 / 2 + z)" would be more traditional*)
   186 *}
   187 ML {*
   188 val (numerator, denominator) = HOLogic.dest_eq eq3;
   189 val eq3' = HOLogic.mk_eq (numerator, fract1); (*A B !*)
   190 term2str eq3';
   191 *}
   192 ML {* (*MANDATORY: otherwise 3 = 0*)
   193 val SOME (eq3'' ,_) = rewrite_set_ @{theory Isac} false norm_Rational eq3';
   194 term2str eq3'';
   195 *}
   196 
   197 subsubsection {*get first koeffizient*}
   198 
   199 ML {*
   200 (*substitude z with the first zeropoint to get A*)
   201 val SOME (eq4_1,_) = rewrite_terms_ @{theory Isac} e_rew_ord e_rls [s_1] eq3'';
   202 term2str eq4_1;
   203 *}
   204 ML {*
   205 val SOME (eq4_2,_) = rewrite_set_ @{theory Isac} false norm_Rational eq4_1;
   206 term2str eq4_2;
   207 *}
   208 ML {*
   209 val fmz = ["equality (3 = 3 * A / (4::real))", "solveFor A","solutions L"];
   210 val (dI',pI',mI') =("Isac", ["univariate","equation"], ["no_met"]);
   211 
   212 *}
   213 ML {*
   214 (*solve the simple linear equilation for A*)
   215 val (p,_,f,nxt,_,pt) = CalcTreeTEST [(fmz, (dI',pI',mI'))];
   216 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   217 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   218 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   219 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   220 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   221 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   222 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   223 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   224 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   225 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   226 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   227 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   228 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   229 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   230 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   231 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   232 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   233 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   234 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   235 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   236 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   237 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   238 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   239 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   240 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   241 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   242 *}
   243 ML {*
   244 val (p,_,f,nxt,_,pt) = me nxt p [] pt; 
   245 f2str f;
   246 *}
   247 
   248 subsubsection {*get second koeffizient*}
   249 
   250 ML {*
   251 (*substitude z with the second zeropoint to get B*)
   252 val SOME (eq4b_1,_) = rewrite_terms_ @{theory Isac} e_rew_ord e_rls [s_2] eq3'';
   253 term2str eq4_1;
   254 *}
   255 
   256 ML {*
   257 val SOME (eq4b_2,_) = rewrite_set_ @{theory Isac} false norm_Rational eq4b_1;
   258 term2str eq4b_2;
   259 *}
   260 
   261 ML {*
   262 (*solve the simple linear equilation for B*)
   263 
   264 val fmz = ["equality (3 = -3 * B / (4::real))", "solveFor B","solutions L"];
   265 val (dI',pI',mI') =("Isac", ["univariate","equation"], ["no_met"]);
   266 
   267 val (p,_,fb,nxt,_,pt) = CalcTreeTEST [(fmz, (dI',pI',mI'))];
   268 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   269 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   270 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   271 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   272 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   273 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   274 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   275 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   276 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   277 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   278 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   279 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   280 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   281 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   282 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   283 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   284 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   285 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   286 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   287 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   288 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   289 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   290 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   291 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   292 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   293 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;
   294 val (p,_,fb,nxt,_,pt) = me nxt p [] pt; 
   295 f2str fb;
   296 *}
   297 
   298 ML {* (*check koeffizients*)
   299 if f2str f = "[A = 4]" then () else error "part.fract. eq4_1";
   300 if f2str fb = "[B = -4]" then () else error "part.fract. eq4_1";
   301 *}
   302 
   303 end
   304