kill(all)$

maxima : Code : Muller method (algo),version 1.0,part 1
peter.vlasschaert@gmail.com, 01/08/2020

see website :
==========
maxima : graphics for Muller method (algo),version 1.0
peter.vlasschaert@gmail.com, 01/08/2020
----------------------------------------------------------------------------
use : plot equation,q1a,q1b,q1c

load packages ,option
=================

load(lrats)$
ratprint:false$

the three starting values : parabolic profile

x0:4.5;
x1:5.5;
x2:5.0;/*ref value */

m:[ 20.625,48,82.875]$

q1a:['parametric, x0, 0+t, [t, 0, m[1]]];
q1b:['parametric, x2, 0+t, [t, 0, m[2]]];
q1c:['parametric, x1, 0+t, [t, 0, m[3]]];

only the parametric lines  of  intersection with nonlinear equation.

wxplot2d([q1a,q1b,q1c],[x,3,6],[y,-40,100],[nticks,200]
,[gnuplot_preamble,"set size ratio 1;"],[box,true],[axes,true],[legend,false])$

general line (beginpoint to endpoint)
=========
x0:4.5,x1:5.5,x2:5

first line :
-----------

 begin point : (x0,m[1])
 mid   point :  (x2,m[2])

rico1:(m[2]-m[1])/(x2-x0);

second line :
-----------------

mid   point  : (x2,m[2])
end   point :  (x1,m[3])

rico2:(m[2]-m[3])/(x2-x1);

How to draw parametric line:(with parameter: t)
====================================
t= 0 , we get first point of the line      , coordinates   (xx1,yy1)
t=1  , we get second point of the line, coordinates  (xx2,yy2)
-----------------------------------------------------------------
draw line : 0 ≤ t ≤ 1 

p1:xx = a1*xx1+a2*xx2;

find a1,a2 ? 
set : a2 = t

p2:lratsubst([a2=t],p1);

find a1 ?

p2a:lratsubst([a1=1-t],p2);

same yy 

p2b:yy=t*yy2+(1-t)*yy1;

qaa:['parametric, rhs(p2a), rhs(p2b), [t, 0, 1]];

mp1:[xx1= x0,xx2 = x2,yy1=m[1],yy2=m[2]];

qa:lratsubst(mp1,qaa);

slope : first line (purple),see below, (see:rico1 ,above)

wxplot2d([q1a,q1b,q1c,qa],[x,3,6],[y,-40,100],[nticks,200]
,[gnuplot_preamble,"set size ratio 1;"],[box,true],[axes,true],[legend,false])$

slope: second line (black), see below (see:rico2 ,above)

mp2:[xx1= x2,xx2 = x1,yy1=m[2],yy2=m[3]];
qb:lratsubst(mp2,qaa);

wxplot2d([q1a,q1b,q1c,qa,qb],[x,3,6],[y,-40,100],[nticks,200]
,[gnuplot_preamble,"set size ratio 1;"],[box,true],[axes,true],[legend,false])$

the way you can construct figure with parametric lines.( only use plot2d)
wx : before , be on the canvas .(wxplot2d(......))

1e)  parabolic equation with (ref :x2 here x[2])
--------------------------------------------------------------

eq[i]:g[a,b,c,xref,i]:=a*(x[i]-xref)^2+b*(x[i]-xref)+c;

for i: 0,1,2

eq[0]:g[a,b,c,x[2],0];
eq[1]:g[a,b,c,x[2],1];
eq[2]:g[a,b,c,x[2],2];

eq[0];

2e)  f()  = nonlinear equation
---------------------------------------
example :f(x) =x^3-13*x-12

s1:f(x[0])=eq[0];
s2:f(x[1])=eq[1];
s3:f(x[2])=eq[2];

a) we need elimate (a & b) , because c = see 's3'

sa1:s1-s3;
sa2:s2-s3;

b) we set x[0],x[1],x[2]

sb1:x[1]-x[0]=h[0];
sb2:x[2]-x[1]=h[1];
sc1: ((f(x[1])-f(x[0]))/lhs(sb1))=d[0];
sc2: ((f(x[2])-f(x[1]))/lhs(sb2))=d[1];

needed for sa1,sa2

sc1a:f(x[1])-f(x[0]) = d[0]*h[0];
sc1b:f(x[2])-f(x[1]) = d[1]*h[1];
sc1c:sc1a+sc1b;

needed for sa1,sa2

sb3:sb2+sb1;

sd1:lhs(sa1)=factorsum(lratsubst([sb3],rhs(sa1)));
sd2:lhs(sa2)=lratsubst([sb2],rhs(sa2));

sdf1:lratsubst([sc1c],sd1)*(-1),factorsum;
sdf2:lratsubst([sc1b],sd2)*(-1);

solve :a,b
--------------

sol1:flatten(solve([sdf1,sdf2],[a,b]));

t1:[rhs(sol1[1])=lhs(sol1[1])];

sol2:lratsubst(t1,sol1[2]);

sol3:lhs(sol2)=factorsum(rhs(sol2));

sol:[sol1[1],sol3];

To draw the graph we need : a,b,c (see : website)
====================================

x[0]:x0;
x[1]:x1;
x[2]:x2;

z1:[ev(sb1,nouns),ev(sb2,nouns)];

ff1:f(x) := x^3-13*x-12;

ps1:lratsubst([x=x0],ff1),numer;
ps2:lratsubst([x=x1],ff1),numer;
ps3:lratsubst([x=x2],ff1),numer;

find  → ' sc1 : value '

sc1aa:ev(sc1,nouns);

find  → ' sc2 : value '

sc1bb:ev(sc2,nouns);

z2:[sc1aa,sc1bb];

z:append(z1,z2);

rem : z can't be used by 'lratsubst' , you get error.

zz:makelist(rhs(z[i])=lhs(z[i]),i,1,length(z));

find a,b :use 'sol' ,above

solf : lratsubst(zz,sol);

solfinal:append(solf,[c=f(x[2])]);

What sign: (discriminant), sign(b)

definition :sign
==========

sign(-2);

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

rationalized :version , see : Appendix : Muller Method

maxima : pratical Muller method (= algo ),version 1.0
peter.vlasschaert@gmail.com,31/07/2020

Δ:lratsubst(solfinal,b^2-4*a*c),numer;

sg:sgn[rhs(solfinal[2])];

xm: lratsubst(solfinal,-2*c/(b+sg*sqrt(Δ))),numer;

the first : root approximation : from allrroots = 4
=================================
→ first iteration to find the root

x[3]=x[2]+xm;