Artifact 06fe165a9feb68b3ffa92a54724c4e1ac23a1baded2935741733fac6790d6a34:
- Executable file
r37/packages/trigint/trigint.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 9534) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/trigint/trigint.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 9534) [annotate] [blame] [check-ins using]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Neil Langmead November 1996 ZIB Berlin % routines to evaluate trigonometric integrals % % main routine is trigint % substitution variable is always u % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% module trigint; create!-package ('(trigint),nil); global '(!*tracetrig); switch tracetrig; off tracetrig; % off by default; algebraic; load_package limits; on usetaylor; load_package misc; % need a knowledge base of the substitutions expr procedure sub_a(exp,var); sub({sin(var)=2*u/(1+u^2), cos(var)=(1-u^2)/(1+u^2)},exp); expr procedure sub_b(exp,var); sub({sin(var)=(u^2-1)/(u^2+1), cos(var)=2*u/(u^2+1)},exp); expr procedure sub_c(exp,var); sub({sin(var)=2*u/(1+u^2), cos(var)=(u^2-1)/(1+u^2)},exp); expr procedure sub_d(exp,var); sub({sin(var)=u/(sqrt(1+u^2)), cos(var)=1/(sqrt(1+u^2))},exp); % applying the substitutions to their integrals expr procedure apply_a(exp,var); begin scalar answer, result; answer:=sub_a(exp,var); answer:=answer* 2/(1+u^2); result:=int(answer,u); result:=sub({u=tan(var/2)},result); result:=result+k(result,var,pi)*floor((var-pi)/(2*pi)); return result; end; expr procedure apply_b(exp,var); begin scalar answer, result; answer:=sub_b(exp,var); answer:=answer*2/(1+u^2); result:=int(answer,u); result:= sub({u=tan(var/2+pi/4)}, result); result:=result+k(result,var,pi/2)*floor((var-pi/2)/(2*pi)); return result; end; expr procedure apply_c(exp, var); begin scalar answer, result; answer:=sub_c(exp,var); answer:= answer*(-2/(1+u^2)); result:=int(answer,u); result:=sub({u=1/(tan(var/2))},result); result:=result+k(result,var,0)*floor(var/(2*pi)); return result; end; expr procedure apply_d(exp, var); begin scalar answer, result; answer:=sub_d(exp,var); answer:=answer*(1/(1+u^2)); result:=int(answer,u); result:= sub({u=tan(var)},result); result:=result+k(result,var,pi/2)*floor((var-pi/2)/pi); return result; end; % some trigonometric substitutions which occur very frequently trig_rules:= { (sec(~x))^2 => 1/(cos(x))^2, %tan(~x) => sin(x)/cos(x), (tan(~x))^2 => (sec(x))^2-1 %2*sin(~x/2)*cos(~x/2) => sin(x) %(cos(~x))^2 => 1-(sin(x))^2 % (sin(~x))^2 => 1-(cos(x))^2 }; %let trig_rules; %for all x let sin(x)^2+cos(x)^2=1, % 2*sin(x/2)*cos(x/2)=sin(x); % procedure to see if integral is returned unevaluated expr procedure unevalp(exp); begin scalar finished, k; k:=0; finished:=0; while ( finished=0 and k<=arglength(exp)) do << if(part(exp,k)=int) then finished:=1 else k:=k+1; >>; if (finished=1) then return t else nil; end; % procedure to evaluate K expr procedure k(exp,var,val); limit!-(exp,var,val)-limit!+(exp,var,val); % two routines to see if we have either unevaluated limits or ints in our % result. If we have, then try a different substitution expr procedure uneval_int(exp); begin scalar temp; if( not freeof(exp,int)) then return temp:=t else temp:=nil; return temp; end; expr procedure uneval_lim(exp); begin scalar temp; if(not freeof(exp,limit!-)) then return t else << if(not freeof(exp,limit!+)) then return t else nil; >>; end; expr procedure fail_a(exp,var); begin scalar temp; temp:=apply_a(exp,var); if(uneval_lim(temp)) then return t else << if(uneval_int(temp)) then return t else return nil; >>; end; expr procedure fail_b(exp,var); begin scalar temp; temp:=apply_b(exp,var); %temp:=temp+k(temp,var,pi/2)*floor((var-pi/2)/(2*pi)); if(uneval_lim(temp)) then return t else << if(uneval_int(temp)) then return t else return nil; >>; end; expr procedure fail_c(exp,var); begin scalar temp; temp:=apply_c(exp,var); %temp:=temp+k(temp,var,0)*floor(var/(pi)); if(uneval_lim(temp)) then return t else << if(uneval_int(temp)) then return t else return nil; >>; end; expr procedure fail_d(exp,var); begin scalar temp; temp:=apply_d(exp,var); %temp:=temp+k(temp,var,pi/2)*floor((var-pi/2)/pi); if(uneval_lim(temp)) then return t else << if(uneval_int(temp)) then return t else return nil; >>; end; expr procedure fail(exp); if(uneval_lim(exp)) then t else << if(uneval_int(exp)) then t else nil; >>; let log(-1) => i*pi; % really important. If further log rules are needed, % take a look at the ratint package, and the module % convert, which contains an extensive list of such % rules expr procedure trigint(exp,var); begin scalar answer, answer_1, answer_2, answer_3, answer_4, result; %off mesgs; % off by default % check for correct input if(freeof(exp,sin) and freeof(exp,cos) and freeof(exp,tan)) then << if(lisp !*tracetrig) then write "expression free of sin, cos tan, proceeding with standard integration"; return int(exp,x); >>; on usetaylor; if freeof(exp,sin(var)) then % we use substitution (a) << answer:=apply_a(exp,var); %answer:=answer+k(answer,var,pi)*floor((var-pi)/(2*pi)); if(fail(answer)) then % system can't evaluate after subs << if(lisp !*tracetrig) then write "system can't integrate after substitution A, trying again"; answer_2:=apply_b(exp,var); if(fail(answer_2)) then << if(lisp !*tracetrig) then write "trying again with substitution B"; answer_3:=apply_c(exp,var); if(fail(answer_3)) then << if(lisp !*tracetrig) then write "and again with substitution C"; answer_4:=apply_d(exp,var); if(fail(answer_4)) then << if(lisp !*tracetrig) then write "failed in all attempts, system cannot integrate"; return answer; >> else return answer_4; >> else return answer_3; >> else return answer_2; >> else return answer; %let trig_rules; %if(unevalp(answer)) then rederr "system cannot integrate after subs" %else nil; >> else % we use substitution b,c or d << if(freeof(exp,cos(var))) then % use substitution b << answer:=apply_b(exp,var); if(fail(answer)) then << if(lisp !*tracetrig) then write "failed with substitution B: system could not integrate after subs, trying A"; answer_2:=apply_a(exp,var); if(fail(answer_2)) then << if(lisp !*tracetrig) then write "failed with A: trying C now"; answer_3:=apply_c(exp,var); if(fail(answer_3)) then << if(lisp !*tracetrig) then write "failed with C: trying D now"; answer_4:=apply_d(exp,var); if(lisp !*tracetrig) then write "trying all possible substitutions"; if(fail(answer_4)) then rederr "system can't integrate after subs" else return answer_4; >> else return answer_3; >> else return answer_2; >> else return answer; >> else << % now describe situations best for (c) and (d) G and R sect 2.504 if(sub({sin(var)=-sin(var),cos(var)=-cos(var)},exp)=exp) then % d is the best sub in this case << if(lisp !*tracetrig) then write "using heuristics: G & R section 2.504 to integrate "; answer:=apply_d(exp,var); if(fail(answer)) then << if(lisp !*tracetrig) then write "subs D failed, trying now with A"; answer_2:=apply_a(exp,var); if(fail(answer_2)) then << if(lisp !*tracetrig) then write "subs B falied, trying with sub C"; answer_3:=apply_b(exp,var); if(fail(answer_3)) then << if(lisp !*tracetrig) then write "sub C falied, trying sub D"; answer_4:=apply_c(exp,var); if(fail(answer_4)) then rederr "can't integrate after subs" else return answer_4; >> else return answer_3; >> else return answer_2; >> else return answer; >> else << % no guidelines, try each substitution in turn, and return the % best possible answer, if there is one answer:=apply_a(exp,var); if(fail(answer)) then << if(lisp !*tracetrig) then write "not using heuristics, attempting subs in order: trying A"; answer_2:=apply_b(exp,var); if(fail(answer_2)) then << if(lisp !*tracetrig) then write "A failed, trying B"; answer_3:=apply_c(exp,var); if(fail(answer_3)) then << answer_4:=apply_d(exp,var); if(fail(answer_4)) then rederr "can't do it" else return answer_4; >> else return answer_3; >> else return answer_2; >> else return answer; >>; >>; % for the else just before the comment on G and R >>; return answer; end; endmodule; end;