maxima :Muller's method, theoretical(example),version 1.0
peter.vlasschaert@gmail.com,22/07/2020

kill(all)$

load packages
===========

load(lrats)$

option 
====

ratprint:false;

Muller Method (nonlinear equations : roots)
===============================

a) setting: Muller method
--------------------------------

m:g[xx] :=a*(xx-x[i])^2+b*(xx-x[i])+c;

m1:[x[i],x[i-1],x[i-2]];

m1a:makelist(part(rhs(m),k,1),k,1,length(m1)-1);
m1b:[part(rhs(m),length(m1))];
m1c:flatten(append(m1a,m1b));

part(m1[2],1);

eq[k]:=f[part(m1[k],1)]=g[m1[k]];

m2:makelist(eq[k],k,1,length(m1));

p1:flatten(solve(m2,m1c));

rem : we need sign of b:pos  to select sign before : sign(b)* sqrt(b^2-4*ac)

b)  when : g[x[i+1]]=0 ,so we can solve  for : h = x[i+1]-x[i]
---------------------------------------------------------------------------

p2:h=x[i+1]-x[i];

p2a:g[x[i+1]]=0;

p2b:lratsubst([rhs(p2)=lhs(p2)],p2a);

which : solution to take , 'see below'

p2c:expand(solve(p2b,h));

p2d::lratsubst([p2],p2c);

p2e:p2d+x[i];

which equation to choose: sgn(b)
========================

sign(-2.0);

sgn[w]:= if neg=sign(w) then
          -1 else +1;

sgn[-2.0];

example:
=======

 ?  find :rr ,solve F[rr]=0
==================

F[rr]:=5/3*cos(40/180*%pi)-5/2*cos(rr)+11/6-cos(40/180*%pi-rr);

evaluation of a function :
=================

F[0.5];

ev(F[0.5],nouns),float;

ev(F[0.5],nouns),numer;

rem :

ev(%pi,float);

ev(%pi,numer);

a) approximations: rr: for find F[rr], muller method , 'three values' 
----------------------------------------------------------------------------------

r1:[30/180*%pi,(30.5)/180*%pi,31/180*%pi];

cos(r1),numer;

b) values
------------
use (m1) :[x[i],x[i-1],x[i-2]]

r2:maplist(":",m1,r1);

r2a:ev(p1,nouns);

c) function values
----------------------
use( m2):[f[i],fi-1],f[i-2]]

r3:makelist(lhs(m2[ks]),ks,1,length(r2));

r3a:makelist(ev(F[r2[ks1]],nouns),ks1,1,length(r2)),numer;

r3b:maplist(":",r3,r3a);

r3c:ev(p1,nouns),numer;

→   math : sign of b

r3c1:sgn[rhs(r3c[2])];

→  math :  sign of b*(sqr(b^2-4*a*c))

why : reason we take sign of b :x[i+1] close to x[i]

d) result
-----------

result :rad ( =radian)
===============

res:lratsubst(r3c,x[i]+(r3c1*sqrt(b^2-4*a*c)-b)/(2*a)),numer;

result:deg(=degree)
===============

res1:res/%pi*180,numer;

→  " check (with :radian) ?  F[res]=0 :"
============================

F[res],numer;