{VERSION 5 0 "IBM INTEL NT" "5.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 223 "G:=(l,u,z)->((1-z)^ l-1)/(ln(1-z)*(1-z)^u):\ngam:=(i,l,u)->coeff(series(G(l,u,z),z,30),z,i ):\nbeta:=(q,j,l,u)->(-1)^(q-u-j)*\n sum('binomial(i, q-u-j)*gam(i,l,u)',\n 'i'=q-u-j..q-u):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 255 "controleer:=proc(methode ::matrix)\nlocal methode2;\n if(methode[1,1] = 0 and methode[2,1] = 0 ) then\n methode2:=linalg[delcols](methode,1..1);\n methode2:=co ntroleer(methode2);\n else\n methode2:=methode;\n end if;\n retu rn evalm(methode2);\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 632 "AdamsBashforth:=proc(R::equation)\n local s,methode,L,q,a,c; \n if (nargs <> 1) then\n ERROR(\"Er is juist \351\351n gegeven \+ verwacht.\");\n else \n a:=lhs(R);\n c:=rhs(R);\n L:=\"\"||a ;\n if (L=\"p\") then\n q:=c;\n elif (L=\"k\") then\n \+ q:=c;\n elif (L=\"q\") then\n q:=c;\n else\n ERROR(\"D e parameter is niet juist: [p,k,q]=constante.\");\n end if;\n me thode:=linalg[matrix](2,q+1);\n for s from 1 to q+1 do\n metho de[1,s]:=0;\n methode[2,s]:=beta(q,s-1,1,1);\n end do;\n me thode[1,q+1]:=1;\n methode[1,q]:=-1;\n methode[2,q+1]:=0;\n r eturn(evalm(controleer(methode)));\n end if;\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 976 "AdamsMoulton:=proc(R::equation)\n \+ local s,methode,L,q,a,c; \n if (nargs <> 1) then\n ERROR(\"Er is juist \351\351n gegeven verwacht.\");\n else\n a:=lhs(R);\n c: =rhs(R);\n L:=\"\"||a;\n if (L=\"p\") then\n q:=c-1;\n e lif (L=\"q\") then\n q:=c;\n else\n q:=c; \n end if;\n if (((L=\"p\") or (L=\"k\")) and (c=0)) then\n ERROR(\"Deze m ethode bestaat niet.\");\n end if;\n if (L=\"k\" and (c=1)) then \n methode[1]:=linalg[matrix](2,2,[[-1,1],[1/2,1/2]]);\n met hode[2]:=linalg[matrix](2,2,[[-1,1],[0,1]]);\n return ([evalm(met hode[1]),evalm(methode[2])]);\n else\n if(q<>0) then\n \+ methode:=linalg[matrix](2,max(2,q+1));\n for s from 1 to q+1 do \n methode[1,s]:=0;\n methode[2,s]:=beta(q,s-1,1,0); \n end do;\n methode[1,q+1]:=1;\n methode[1,q]:=- 1;\n else\n methode:=linalg[matrix](2,2,[[-1,1],[0,1]]);\n end if;\n return(evalm(controleer(methode)));\n end if; \n end if;\nend proc:" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 766 "Nystrom:=proc(R::equation)\n local s,methode,L,q,a, c; \n if (nargs <> 1) then\n ERROR(\"Er is juist \351\351n gegeve n verwacht.\");\n else\n a:=lhs(R);\n c:=rhs(R);\n L:=\"\"|| a;\n if (L=\"p\") then\n q:=c;\n elif (L=\"q\") then\n \+ q:=c;\n else\n q:=c; \n end if;\n if (((L=\"k\") and (c =1)) or ((L=\"p\") and (c=1))) then\n ERROR(\"Deze methode bestaa t niet.\");\n end if;\n if (q = 1) then\n methode:=linalg[m atrix](2,3,[[-1,0,1],[0,2,0]]);\n else\n methode:=linalg[matri x](2,q+1);\n for s from 1 to q+1 do\n methode[1,s]:=0;\n \+ methode[2,s]:=beta(q,s-1,2,1);\n end do;\n methode[1,q +1]:=1;\n methode[1,q-1]:=-1;\n methode[2,q+1]:=0;\n end \+ if;\n return(evalm(controleer(methode)));\n end if;\nend proc:" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1199 "MilneSimpson:=proc(R::equ ation)\n local s,methode,L,q,a,c; \n if (nargs <> 1) then\n ERRO R(\"Er is juist \351\351n gegeven verwacht.\");\n else\n a:=lhs(R) ;\n c:=rhs(R);\n L:=\"\"||a;\n if (L=\"p\") then\n q:=c- 1;\n elif (L=\"q\") then\n q:=c;\n else\n q:=c; \n \+ end if;\n if (((L=\"k\") and ((c=0) or (c=1))) or ((L=\"p\") and (( c=0) or (c=3)))) then\n ERROR(\"Deze methode bestaat niet.\");\n \+ end if;\n if (L=\"k\" and (c=2)) then\n methode[1]:=linalg[ matrix](2,3,[[-1, 0, 1], [1/3, 4/3, 1/3]]);\n methode[2]:=linalg[ matrix](2,3,[[-1, 0, 1], [0, 2, 0]]);\n methode[3]:=linalg[matrix ](2,3,[[-1, 0, 1], [0, 0, 2]]);\n return([evalm(methode[1]),evalm (methode[2]),evalm(methode[3])]);\n else\n if (q = 0) then\n \+ methode:=linalg[matrix](2,3,[[-1,0,1],[0,0,2]]);\n elif (q \+ = 1) then\n methode:=linalg[matrix](2,3,[[-1,0,1],[0,2,0]]);\n \+ else\n methode:=linalg[matrix](2,max(3,q+1));\n for s from 1 to q+1 do\n methode[1,s]:=0;\n methode[2,s ]:=beta(q,s-1,2,0);\n end do;\n methode[1,q+1]:=1;\n \+ methode[1,q-1]:=-1;\n end if;\n return(evalm(controleer( methode)));\n end if;\n end if;\nend proc:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 626 "BDF:=proc(R::equation)\nlocal j, delt, beta, \+ methode,L,q,a,c;\nif(nargs <> 1) then\n ERROR(\"Er is juist \351\351n gegeven verwacht.\"); \nelse\n a:=lhs(R);\n c:=rhs(R);\n L:=\"\"| |a;\n if (L=\"p\") then\n q:=c;\n elif (L=\"q\") then\n q:=c; \n else\n q:=c; \n end if;\n methode:=linalg[matrix](2,q+1);\n \+ delt:=i->subs(t=0, diff((-1)^i*expand(binomial(-t,i)),t));\n beta:=1/ sum('delt(i)','i'=0..q); \n for j from 0 to q-1 do\n methode[1,j+1 ]:=sum('(-1)^(q-j)*binomial(i,q-j)*delt(i)','i'=q-j..q)*beta;\n met hode[2,j+1]:=0;\n end do;\n methode[1,q+1]:=1:\n methode[2,q+1]:=be ta;\n return(evalm(controleer(methode)));\nend if;\nend proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 910 "Maximaal:=proc(k::posint)\n local alpha,beta,methode,u,rho,i,l,eqn,eqns,vars,sols,y,f;\nif(nargs<> 1) then\n ERROR(\"Er is juist \351\351n gegeven verwacht.\");\nelif ( irem(k,2)<>0) then\n ERROR(\"Deze procedure werkt enkel met een even \+ aantal stappen.\");\nelse\n alpha:=array(0..k);\n beta:=array(0..k); \n u:=array(0..k/2-1);\n rho:=expand((r^2-1)*\n product((r-u[ i])^2+1-(u[i])^2,'i'=1..(k-2)/2));\n for i from 0 to k do alpha[i]:=c oeff(rho,r,i) od;\n eqn:=array(0..k+3);\n y := x-> x^l;\n f := x-> \+ x^(l-1)*l;\n for l from 0 to k+3 do \n eqn[l]:=sum('alpha[i]*y(x0+ i*h)\n -h*beta[i]*f(x0+i*h)','i'=0..k);\n end do;\n eqn s:=\{seq(eqn[l],l=0..k+2)\};\n vars:=\{seq(beta[l],l=0..k)\};\n sols :=solve(eqns,vars);\n assign(sols);\n methode:=linalg[matrix](2,k+1) ;\n for i from 1 to k+1 do\n methode[1,i]:=alpha[i-1];\n method e[2,i]:=beta[i-1];\n end do;\n return(evalm(controleer(methode)));\n end if;\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 169 "uits chrijven:=proc(methode::matrix)\nlocal k;\nk:=linalg[coldim](methode); \nreturn(sum(methode[1,k+1-i]*y[n-i+2],i=1..k)=h*sum(methode[2,k+1-i]* f[n-i+2],i=1..k));\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 493 "orde:=proc(methode::matrix)\n local lmmorde,expr,continu,k;\n \+ k:=linalg[coldim](methode)-1;\n lmmorde:=0:\n continu:=true:\n exp r:=sum(methode[1,i+1]*i^(lmmorde)/lmmorde!,'i'=0..k):\n if(expr=0) th en \n while continu do\n lmmorde:=lmmorde+1;\n expr:=sum( methode[1,i+1]*(i+1)^lmmorde/lmmorde! -\n methode[2,i+1 ]*(i+1)^(lmmorde-1)/(lmmorde-1)!,'i'=0..k);\n if (expr<>0) then \+ \n continu:=false; \n end if;\n end do;\n end if;\n r eturn(lmmorde-1);\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 283 "foutconstante:=proc(methode::matrix)\n local lmmorde,expr,k;\n k:=linalg[coldim](methode)-1;\n lmmorde:=orde(methode):\n expr:=su m(methode[1,i+1]*(i+1)^(lmmorde+1)/(lmmorde+1)! - \+ methode[2,i+1]*(i+1)^lmmorde/lmmorde!,'i'=0..k);\n return(e xpr);\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 152 "PLAF:= proc(methode::matrix)\n local lmmorde;\n lmmorde:=orde(methode):\n \+ return(foutconstante(methode)*h^(lmmorde+1)*(D@@(lmmorde+1))(y)(xi)); \nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 239 "stabvt:=pro c(methode::matrix)\n local veelterm,q;\n q:=linalg[coldim](methode)- 1;\n veelterm:=sum('methode[1,j+1]*r^j','j'=0..q)-h*\n s um('methode[2,j+1]*r^j','j'=0..q);\n return(sort(collect(veelterm,r), [r,h],plex));\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 201 "randkromme:=proc(methode::matrix)\n local veelterm, hbound;\n v eelterm:=stabvt(methode);\n hbound:=solve(veelterm,h);\n print(theta = 0..2*Pi);\nRETURN(subs(r=cos(theta)+I*sin(theta),hbound));\nend pro c:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1993 "plotten:=proc(veelt erm::polynom,GRENSWAARDE::realcons)\nlocal GRENS, hbound, loopbound, A ANTALPUNTEN, nietpunt,plot_nietarcering, bound,i,j,k,l,n,plot_gebied,p lot_arcering,moduli,maxmodulus,punt,wortels,z,zmin,zmax;\n if(nargs = 2) then \n GRENS:= abs(GRENSWAARDE);\n else\n GRENS:= 0.05;\n \+ end if;\n hbound:=solve(veelterm,h);\n assume(x,real):\n bound:=si mplify(subs(r=cos(x)+I*sin(x),hbound));\n plot_gebied:=plot([Re(bound ),Im(bound),x=0..2*Pi],style =line,color=blue);\n k:=1;\n AANTALPUNT EN:=30;\n punt:=[];\n while (nops(punt)=0 and GRENS > 0.001) do\n \+ nietpunt:=[];\n l:=1;\n for i from 0 to AANTALPUNTEN do\n \+ loopbound:=subs(x=evalf(i*Pi/AANTALPUNTEN),bound);\n z:= evalf(Re (loopbound)+I*Im(loopbound));\n zmin:= z-GRENS;\n zmax:= z+G RENS;\n wortels:=fsolve(subs(h=zmin,veelterm)=0,r,complex);\n \+ moduli:=map(abs,[wortels]); \n maxmodulus:=sort(moduli)[-1];\n \+ if (maxmodulus < 1) then\n punt:=[op(punt),evalf([Re(zmin) ,Im(zmin)]),evalf([Re(zmin),-Im(zmin)])];\n k:=k+2;\n else \n nietpunt:=[op(nietpunt),evalf([Re(zmin),Im(zmin)]),evalf([R e(zmin),-Im(zmin)])];\n l:=l+2;\n wortels:=fsolve(subs(h =zmax,veelterm)=0,r,complex);\n moduli:=map(abs,[wortels]); \n \+ maxmodulus:=sort(moduli)[-1];\n if (maxmodulus < 1) then \n punt:=[op(punt),evalf([Re(zmax),Im(zmax)]),evalf([Re(zmax) ,-Im(zmax)])];\n k:=k+2;\n else\n nietpunt:=[ op(nietpunt),evalf([Re(zmax),Im(zmax)]),evalf([Re(zmax),-Im(zmax)])]; \n l:=l+2;\n end if; \n end if;\n end do;\n \+ GRENS:=GRENS/10;\n end do;\n if(nops(punt)<>0) then\n plot_arcer ing:=plots[listplot]([seq(punt[n],n=1..k-1)],style=point,color=green); \n end if;\n plot_nietarcering:=plots[listplot]([seq(nietpunt[n],n=1 ..l-1)],style=point,color=red);\n if(nops(punt)=0) then\n plots[di splay](\{plot_gebied,plot_nietarcering\});\n else \n plots[display ](\{plot_gebied,plot_arcering,plot_nietarcering\});\n end if;\nend pr oc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 636 "stabgebied:=proc(me thode::matrix,GRENSWAARDE::realcons)\n local veelterm,GRENS,plot_niet arcering,plot_gebied,plot_arcering;\n if (linalg[equal](methode,linal g[matrix]([[-1,1],[1/2,1/2]]))) then\n plot_gebied:=plot([0,t,t=-5. .5],style =line,color=blue);\n plot_arcering:=plot([-1,t,t=-5..5],s tyle=point,color=green);\n plot_nietarcering:=plot([1,t,t=-5..5],st yle=point,color=red);\n plots[display](\{plot_gebied,plot_arcering, plot_nietarcering\});\n else\n if(nargs = 2) then \n GRENS:= \+ abs(GRENSWAARDE);\n else\n GRENS:= 0.05;\n end if;\n vee lterm:=stabvt(methode);\n plotten(veelterm,GRENS);\n end if;\nend \+ proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1236 "stabinterval:=pr oc(methode::matrix) \nlocal goed,vt,grens,step,j,ok,opl,tol ;\ngoed:=t rue;\nfor j from 7 to 15 do\n if (linalg[equal](methode,BDF(k=j))) th en\n ERROR(\"Deze methode heeft geen stabiliteitsinterval dat begin t bij 0 (<).\");\n end if;\nend do;\nfor j from 1 to 6 do\n if (lina lg[equal](methode,BDF(k=j))) then\n goed:=false\n end if;\nend do; \nif goed then\n if (linalg[equal](methode,linalg[matrix]([[-1,1],[0, 1]]))) then\n printf(`stabiliteitsinterval : `);\n print([-infin ity,0],U,[2,infinity]);\n elif (linalg[equal](methode,linalg[matrix]( [[-1,1],[1/2,1/2]]))) then\n printf(`stabiliteitsinterval : `);\n \+ print([-infinity,0]);\n else\n vt:=stabvt(methode);\n tol:=1/ 10.^12;\n step:=0.2;\n j:=0;\n ok:=true;\n grens:=0;\n \+ while step>tol do\n while ok do\n j:=j+1; \n opl:=s olve(subs(h=-(grens+j*step),vt)=0,r);\n if sort(map(abs,[opl])) [-1]>1 then \n ok:=false; \n end if;\n end do;\n \+ grens:=grens+(j-1)*step;\n j:=0;\n ok:=true;\n ste p:=step/10.;\n end do; \n printf(`stabiliteitsinterval :`);\n \+ return([-grens,0]);\n end if;\nelse\n ERROR(\"De methode is een BDF -methode (k<7) en heeft een oneindig stabiliteitsgebied.\");\nend if; \nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1317 "nulstabiel :=proc(methode::matrix)\nlocal q,s,rho_, wortels,wortel,moduli,maxmodu lus,i,imax,goed;\nrho_:=sort(sum('methode[1,j]*r^(j-1)','j'=1..linalg[ coldim](methode)),r,plex);\nprint(rho(r)=rho_);\nq:=degree(rho_,r);\nw ortel:=[fsolve(rho_=0,r,complex)]:\nprint(wortels= wortel);\nwortels:= wortel; \nmoduli:=map(abs,wortels):\nmaxmodulus:=sort(moduli)[-1]:\nim ax:=1:\nwhile abs(wortels[imax])<>maxmodulus do \n imax:=imax+1: \nen d do:\ni:=0:\ngoed:=true:\nif (maxmodulus>1) then\n printf(`niet nuls tabiel : `);\n if type(wortels[imax],realcons) then\n printf(`|%6f | > 1`,wortels[imax])\n else\n if Im(wortels[imax])>0 then \n \+ s:=`+`\n else \n s:=`-` \n end if;\n printf(`|%6f %s % 6f I| > 1)`,\n Re(wortels[imax]),s,abs(Im(wortels[imax])))\n end i f;\nelse\n while goed and (i0 then \n s:=`+` \+ \n else \n s:=`-` \n end if;\n printf(`r=%6f %s \+ %6f I`);\n end if;\n printf(`is dubbele wortel met modulus 1`); \n end if;\nend if;\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2633 "helpLMM:=proc()\nprint(`Alle methodes worden weergegeven en \+ gebruikt als 2-dimensionale matrix. Als de methode de volgende is:\\n \+ `);\nprint(sum(alpha[i]*y[n+i],i=0..k)=h*sum(beta[i]*f[n+i],i=0..k),`` );\nprint(`dan wordt de methode voorgesteld door:\\n `);\nprint(evalm( matrix(2,4,[[alpha[0],alpha[1],~,alpha[k]],[beta[0],beta[1],~,beta[k]] ])));\nprint(`Besproken procedures: \\`AdamsBashforth\\`, \\`AdamsMoul ton\\`, \\`MilneSimpson\\`, \\`Nystrom\\`, \\`BDF\\`, \\`Maximaal\\`, \+ \\`uitschrijven\\`, \\`orde\\`, \\`PLAF\\`, \\`stabvt\\`, \\`grens\\`, \\`stabgebied\\`, \\`stabinterval\\`, \\`nulstabiel\\`.\\n`);\nprint( `Bij de constructieprocedures (\\`AdamsBashforth\\`, \\`AdamsMoulton\\ `, \\`MilneSimpson\\`, \\`Nystrom\\` en \\`BDF\\`) is de enige paramet er een vergelijking van de vorm:\\n`);\nprint(p= ~ ,`(Geef de orde van de methode), `);\nprint(k= ~ ,`(Geef het aantal stap pen van de methode), `);\nprint( q= ~ ,`(Geef de constructieparameter (zie cursus)).`);\nprint(`Aan \\`Maximaal\\` \+ geef je \351\351n parameter mee: een positief even getal -> het stapge tal. Als output krijg je telkens de gevraagde methode (in matrixvorm). \\`uitschrijven\\` geeft de gewone vorm van de methode terug. \\`orde \\` en \\`PLAF\\` spreken voor zich,\\`stabvt\\` geeft de stabiliteits veelterm. \\`grens\\` bepaalt een analytische uitdrukking voor de gren s van het stabiliteitsgebied en \\`stabgebied\\` geeft een plot van da t gebied, waarbij de groene punten in het stabiliteitsgebied liggen en de rode erbuiten. Als input geef je de methode. Bij deze laatste is e r nog een optionele parameter, nl. de afstand tot de randkromme voor p unten die moeten gecontroleerd worden. \\`randkromme\\` geeft een anal ytische uitdrukking (in poolco\366rdinaten) van de randkromme van het \+ stabiliteitsgebied, \\`stabinterval\\` geeft het stabiliteitsinterval \+ en \\`nulstabiel\\` controleert of een methode nulstabiel is. Als inpu t geef je de methode als matrix.\\n`);\nprintf(` AdamsBashforth:=p roc(R::equation)\\n`);\nprintf(` AdamsMoulton:=proc(R::equation)\\ n`);\nprintf(` MilneSimpson:=proc(R::equation)\\n`);\nprintf(` \+ Nystrom:=proc(R::equation)\\n`);\nprintf(` BDF:=proc(R::equation) \\n`);\nprintf(` Maximaal:=proc(k::posint)\\n`);\nprintf(` ord e:=proc(methode::matrix)\\n`);\nprintf(` uitschrijven:=proc(method e::matrix)\\n`);\nprintf(` PLAF:=proc(methode::matrix)\\n`);\nprin tf(` stabvt:=proc(methode::matrix)\\n`);\nprintf(` grens:=proc (methode::matrix)\\n`);\nprintf(` stabgebied:=proc(methode::matrix ,GRENSWAARDE::realcons)\\n`);\nprintf(` stabinterval:=proc(methode ::matrix)\\n`);\nprintf(` nulstabiel:=proc(methode::matrix)\\n`); \nreturn(NULL);\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 482 "exact:=proc (probleem,x)\n local z,ODE,eqns,conds,opl;\n z:=pro c (x) local i; \n RETURN(seq(u[i](x),i=1..linalg[coldim](evalm(pro bleem['matrix'])))); \n end:\n ODE:=convert(evalm(map(u->diff(u,v ),\n [z(v)])-linalg[multiply](evalm(probleem['matrix']),[z(v)])),l istlist):\n eqns:=map(u->u=0,ODE): \n conds:=map(u->u=0,convert(eval m([z(probleem['x0'])]-probleem['y0']),listlist)):\n opl:=dsolve(\{op( eqns)\} union \{op(conds)\},\{z(v)\}):\n RETURN (subs(v=x,subs(opl,[z (v)]))):\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 269 "genereer _beginwaarden:=proc(methode::matrix,probleem,h_)\nlocal i,k,y,opl;\ny[ 1]:=probleem['y0'];\nk:=linalg[coldim](methode)-1;\nopl:=exact(problee m,x);\nfor i from 2 to k do\n y[i]:=eval(subs(\{x=probleem['x0']+(i-1 )*h_\},opl));\nend do;\nreturn([seq(y[i],i=1..k)]);\nend proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 842 "take_step:=proc(n::posint,m ethode::matrix,probleem,h_,y)\nlocal expr,k,i,f,onbek,sol;\nk:=linalg[ coldim](methode)-1;\nf:=(y)->linalg[multiply](probleem['matrix'],y);\n expr:=array(sparse,1..linalg[coldim](probleem['matrix']));\nfor i from 1 to k do\n expr:=linalg[matadd](expr,linalg[matadd](methode[1,i]*y[ n-k+i-1],\n h_*methode[2,i]*f(y[n-k+i-1]),-1 ,1));\nend do;\nif(methode[2,k+1]<>0) then\n onbek:=array(1..linalg[c oldim](probleem['matrix']),1..1);\n sol:=solve(\{seq(onbek[i,1]=(lina lg[matadd](expr,h_*methode[2,k+1]*f(onbek)))[i,1],\n i=1..linalg[ coldim](probleem['matrix']))\},\n \{seq(onbek[i,1],i=1..linalg[co ldim](probleem['matrix']))\});\n assign(sol);\n return([seq(y[i],i=1 ..n-1),[seq(onbek[i,1],i=1..linalg[coldim](probleem['matrix']))]]);\ne lse\n return([seq(y[i],i=1..n-1),evalm(expr)]);\nend if;\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 337 "take_step2:=proc(n::posi nt,methode::matrix,probleem,h_)\nlocal k,y,j;\nk:=linalg[coldim](metho de)-1;\nif (n <= k) then\n ERROR(`uw stap is kleiner dan het aantal b eginwaarden`);\nelse\n y:=genereer_beginwaarden(methode,probleem,h_); \n for j from k+1 to n do\n y:=take_step(j,methode,probleem,h_,y); \n end do;\n return(y);\nend if;\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 313 "LAF:=proc(methode::matrix,probleem,t::posint)\n local y,k,l;\nif (nargs=3) then \n l:=t;\nelse\n l:=orde(methode)+3; \nend if;\nk:=linalg[coldim](methode)-1;\ny:=take_step2(k+1,methode,pr obleem,'h');\nreturn(map(x->series(x,h=0,l),\n evalm(array([op(subs(x =probleem['x0']+(k)*'h',exact(probleem,x)))])-y[k+1])));\nend proc:" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 417 "num:=proc(np::posint,h::r ealcons,methode::matrix,probleem)\nlocal y,n,k,i,plot;\ny:=take_step2( np,methode,probleem,h);\nfor i from 1 to linalg[coldim](probleem['matr ix']) do\nplot[i]:=plots[listplot]([seq([j*h,y[j][i]],j=1..np)],style= point,color=COLOR(RGB, rand()/10^12, rand()/10^12, rand()/10^12)); \+ \nend do;\nplots[display](\{seq(plot[i],i=1. .linalg[coldim](probleem['matrix']))\});\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 479 "GAF:=proc(np::posint,h::realcons,methode ::matrix,probleem)\nlocal y, opl,i,plot;\ny:=take_step2(np,methode,pro bleem,h);\nopl:=exact(probleem,x);\nfor i from 1 to linalg[coldim](pro bleem['matrix']) do\nplot[i]:=plots[listplot]([seq([j*h,subs(x=problee m['x0']+j*h,opl[i])-y[j+1][i]],j=1..np-1)],style=point,color=COLOR(RGB , rand()/10^12, rand()/10^12, rand()/10^12)); \+ end do;\nplots[display](\{seq(plot[i],i=1..linalg[coldim](problee m['matrix']))\});\nend proc:" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2717 "helpOP_PROBLEEM_LMM:=proc();\nprint(`Dezelfde af spraken i.v.m. de input van methodes blijft gelden (typ \\`helpLMM( ) \\` voor meer uitleg).\\n`);\nprint(`Besproken procedures: \\`exact\\` , \\`genereer_beginwaarden\\`, \\`take_step\\`, \\`take_step2\\`, \\`n um\\`, \\`LAF\\` en \\`GAF\\`.\\n`);\nprint(`De procedure \\`exact\\` \+ berekent de exacte oplossing van een lineair probleem. Als het problee m hetvolgende is:\\n`);\nprint(diff(y(x),x) = A*y(x));\nprint(y(x[0]) \+ = y[0]);\nprint(` dan wordt dit voorgesteld door een lijst met drie ar gumenten, nl. \\`matrix\\`, \\`x0\\` en \\`y0\\`.\\n`);\nprintf(` \+ probleem['matrix']:=matrix(m,n,[[~],~,[~]]);\\n`);\nprintf(` probl eem['x0']:=0;\\n`);\nprintf(` probleem['y0']:=[~];\\n`);\nprint(`O pgelet: zowel de quotes rond de namen \\`matrix\\`, \\`x0\\` en \\`y0 \\`, als de namen zelf mogen niet veranderd worden. Om \\`exact\\` dan toe te passen, doe je hetvolgende:\\n`);\nprintf(` exact(probleem ,x);\\n`);\nprint(`Deze \\`x\\` is de variabele waarin men \\`y\\` zal uitdrukken. Belangrijk om weten is dat de volgende procedures enkel w erken met expliciete methodes. De procedure \\`genereer_beginwaarden\\ ` berekent het aantal (exacte) beginwaarden dat nodig is om de methode te kunnen toepassen. Als input geef je de methode, het probleem en de gewenste \\`h\\`-waarde. Men kan ook het symbool \\`h\\` ingeven. De \+ procedure \\`take_step\\` berekent y[x0 + (n-1) h] a.d.h.v een lijst [ y[x0], ... ,y[x0 + (n-2) h]]. Als input geef je dus de \\`n\\`, de met hode, het probleem, de \\`h\\`-waarde en voorgenoemde lijst. De proced ure \\`take_step2\\` doet hetzelfde, maar begint van nul. Het genereer t zelf de beginwaarden en berekent de gevraagde waarden tot het gewens te punt: \\`x0 + (n-1) h\\`. Men moet dus enkel de \\`n\\`, de methode , het probleem en de \\`h\\`-waarde ingeven. \\`LAF\\` berekent een re eksontwikkeling van de lokale afknottingsfout in \\`h\\`. Men geeft de methode en het probleem mee, en optioneel kan nog een derde parameter gebruikt worden om zelf de orde van de ontwikkeling te geven. \\`num \\` geeft een plot van de numerieke oplossing en \\`GAF\\` een plot va n de globale afknottingsfout. Voor beide heeft men dezelfde parameters : het aantal te berekenen punten, de \\`h\\`-waarde, de methode en het probleem.\\n`);\nprintf(` exact:=proc (probleem,x)\\n`);\nprintf( ` genereer_beginwaarden:=proc(methode::matrix,probleem,h)\\n`);\np rintf(` take_step:=proc(n::posint,methode::matrix,probleem,h,y)\\n `);\nprintf(` take_step2:=proc(n::posint,methode::matrix,probleem, h)\\n`);\nprintf(` LAF:=proc(methode::matrix,probleem,t::posint)\\ n`);\nprintf(` num:=proc(n::posint,h::realcons,methode::matrix,pro bleem)\\n`);\nprintf(` GAF:=proc(n::posint,h::realcons,methode::ma trix,probleem)\\n`);\nreturn(NULL);\nend proc:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 93 "rho:=proc(methode::polynom)\nlocal PI;\nPI:=co llect(methode,h);\nRETURN(tcoeff(PI,h));\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "sigma:=proc(methode::polynom)\nlocal PI; \nPI:=collect(methode,h);\nRETURN(-lcoeff(PI,h));\nend proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 244 "zetstapgelijk:=proc(P::matr ix,C::matrix)\nlocal PI_P,PI_C,p,c;\nPI_P:=stabvt(P);\nPI_C:=stabvt(C) ;\np:=degree(PI_P,r);\nc:=degree(PI_C,r);\nif (p>c) then\n PI_C:=r^(p -c)*PI_C;\nelif (c>p) then\n PI_P:=r^(c-p)*PI_P;\nend if;\nRETURN(PI_ P,PI_C);\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 325 "sim ple:=proc(PI::polynom(anything,r))\nlocal R,i,expr;\nR:=op(factor(PI)) ;\nif (type(factor(PI),`*`)) then\n expr:=1;\n for i from 1 to nops( [R]) do\n if (degree(R[i],r)>0 and rem(R[i],r,r)<>0) then\n ex pr:=expr*R[i]\n end if; \n end do;\nelse expr:=PI;\nend if;\nRETU RN(sort(collect(expr,r,factor),[r,h],plex));\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1599 "stabvtPC:=proc(PCpaar)\nlocal M_m u,expr,beta_k,H,PI,W,P_,C_,mu_,t_,L_;\nP_:=PCpaar['P'];\nC_:=PCpaar['C '];\nmu_:=PCpaar['mu'];\nt_:=PCpaar['t'];\nL_:=PCpaar['L'];\nbeta_k:=C _[2,linalg[coldim](C_)];\nif ((L_<>\"Lin\") and (L_<>\"n\") and (L_<> \"Lout\")) then\n ERROR(\"De laatste parameter in PCpaar kan enkel \\ \"n\\\", \\\"Lin\\\" of \\\"Lout\\\" zijn.\");\nend if; \nif ((mu_=0) \+ and (L_<>\"n\")) then\n ERROR(\"Men kan enkel lokale extrapolatie toe passen als men minstens \351\351nmaal de corrector t oepast\");\nend if;\nif (mu_=0) then\n RETURN(stabvt(P_));\nend if;\n H:=beta_k*h;\nM_mu:=x->x^mu_*(1-x)/(1-x^mu_);\nPI:=zetstapgelijk(P_,C_ );\nif (t_<>0 and t_<>1) then\n ERROR(\" 't' kan enkel 0 of 1 zijn.\" )\nend if;\nif (L_=\"n\") then\n if (t_=0) then\n expr:=PI[2]+M_mu (H)*PI[1];\n else\n expr:=beta_k*r^max(linalg[coldim](P_)-1,linalg [coldim](C_)-1)*PI[2]+\n simplify(M_mu(H))*(rho(PI[1])*sigma( PI[2])-rho(PI[2])*sigma(PI[1]));\n end if;\nelse\n W:=foutconstante( C_)/(foutconstante(P_)-foutconstante(C_));\n if (L_=\"Lin\") then\n \+ if (t_=0) then\n expr:=(1+W)*PI[2]+(M_mu(H+W*H)-W)*PI[1] \n \+ else\n expr:=beta_k*r^max(linalg[coldim](P_)-1,linalg[coldim](C_) -1)*((1+W)*PI[2]-W*PI[1])\n +M_mu(H+W*H)*(rho(PI[1])*sigma( PI[2])-rho(PI[2])*sigma(PI[1]));\n end if; \n elif (L_=\"Lout\" ) then\n if (t_=0) then\n expr:=(1+W)*PI[2]+(M_mu(H)+W*H-W)*PI [1]\n else\n expr:=beta_k*r^max(linalg[coldim](P_)-1,linalg[co ldim](C_)-1)*((1+W)*PI[2]-W*PI[1])\n +(M_mu(H)+W*H)*(rho(PI [1])*sigma(PI[2])-rho(PI[2])*sigma(PI[1]));\n end if; \n end if ;\nend if;\nRETURN(simple(expr));\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1022 "ordePC:=proc(PCpaar)\n local p_,p,ORDE,ORDE_n ,laf,P_,C_,mu_,t_,L_,PC,probl;\n P_:=PCpaar['P'];\n C_:=PCpaar['C']; \n mu_:=PCpaar['mu'];\n t_:=PCpaar['t'];\n L_:=PCpaar['L'];\n if ( (L_<>\"Lin\") and (L_<>\"n\") and (L_<>\"Lout\")) then\n ERROR(\"De laatste parameter in PCpaar kan enkel \\\"n\\\", \\\"Lin\\\" of \\\"L out\\\" zijn.\");\n end if; \n if ((mu_=0) and (L_<>\"n\")) then\n \+ ERROR(\"Men kan enkel lokale extrapolatie toepassen als men minstens \351\351nmaal de corrector\n toepast\");\n end if;\n p:=ord e(C_);\n p_:=orde(P_);\n if ((p_>=p) or ((p_p-p_))) the n\n ORDE_n:=p;\n elif ((p_ " 0 "" {MPLTEXT 1 0 570 "randkrommePC:=proc(PCpaar) \nlocal i, veelterm, hbound, L_;\n L_:=PCpaar['L'];\n if ((L_<>\"Lin \") and (L_<>\"n\") and (L_<>\"Lout\")) then\n ERROR(\"De laatste p arameter kan enkel \\\"n\\\", \\\"Lin\\\" of \\\"Lout\\\" zijn.\");\n \+ end if; \n if ((mu_=0) and (L_<>\"n\")) then\n ERROR(\"Men kan en kel lokale extrapolatie toepassen als men minstens \351\351nmaal de co rrector toepast\");\n end if;\n veelter m:=stabvtPC(PCpaar);\n hbound:=solve(veelterm,h);\n print(theta = 0. .2*Pi);\n RETURN(seq(subs(r=cos(theta)+I*sin(theta),hbound[i]),i=1..n ops([hbound])));\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2863 "stabgebiedPC:=proc(PCpaar,GRENSWAARDE::realcons)\nlocal veelterm ,GRENS, hbound, loopbound, AANTALPUNTEN, nietpunt,plot_nietarcering, b ound,i,j,k,l,n,plot_gebied,plot_arcering,moduli,maxmodulus,punt,wortel s,z,zmin,zmax,L_;\n L_:=PCpaar['L'];\n if ((L_<>\"Lin\") and (L_<>\" n\") and (L_<>\"Lout\")) then\n ERROR(\"De laatste parameter kan en kel \\\"n\\\", \\\"Lin\\\" of \\\"Lout\\\" zijn.\");\n end if; \n if ((mu_=0) and (L_<>\"n\")) then\n ERROR(\"Men kan enkel lokale extr apolatie toepassen als men minstens \351\351nmaal de corrector \+ toepast\");\n end if;\n veelterm:=s tabvtPC(PCpaar);\n if(nargs = 2) then \n GRENS:= abs(GRENSWAARDE); \n else\n GRENS:= 0.5;\n end if;\n if (degree(veelterm,h)=1) the n\n plotten(veelterm,GRENS);\n else\n hbound:=solve(veelterm,h) ;\n assume(x,real):\n for i from 1 to nops([hbound]) do\n b ound[i]:=simplify(subs(r=cos(x)+I*sin(x),hbound[i]));\n end do;\n \+ plot_gebied:=plot(\{seq([Re(bound[i]),Im(bound[i]),x=0..2*Pi],i=1..n ops([hbound]))\},\n style=point,symbol=cross,color=bl ue);\n k:=1;\n AANTALPUNTEN:=6;\n punt:=[];\n nietpunt:=[] ;\n while (k=1 and GRENS > 0.001) do\n l:=1;\n for i from 0 to AANTALPUNTEN do\n for j from 1 to nops([hbound]) do\n \+ loopbound:=subs(x=evalf(i*Pi/AANTALPUNTEN),bound[j]);\n \+ z:= evalf(Re(loopbound)+I*Im(loopbound));\n zmin:= z-GRENS; \n zmax:= z+GRENS;\n wortels:=fsolve(evalf(subs(h=zm in,veelterm))=0,r,complex);\n moduli:=map(abs,[wortels]); \n \+ maxmodulus:=sort(moduli)[-1];\n if (maxmodulus < 1) \+ then\n punt:=[op(punt),evalf([Re(zmin),Im(zmin)]),evalf([Re (zmin),-Im(zmin)])];\n k:=k+2;\n else \n \+ nietpunt:=[op(nietpunt),evalf([Re(zmin),Im(zmin)]),evalf([Re(zmin),- Im(zmin)])];\n l:=l+2;\n wortels:=fsolve(evalf(s ubs(h=zmax,veelterm))=0,r,complex);\n moduli:=map(abs,[wort els]); \n maxmodulus:=sort(moduli)[-1]; \n if ( maxmodulus < 1) then\n punt:=[op(punt),evalf([Re(zmax),Im (zmax)]),evalf([Re(zmax),-Im(zmax)])];\n k:=k+2;\n \+ else\n nietpunt:=[op(nietpunt),\n \+ evalf([Re(zmax),Im(zmax)]),evalf([Re(zmax),-Im(zmax)])];\n \+ l:=l+2;\n end if; \n end if;\n end do;\n end do;\n GRENS:=GRENS/10;\n end do;\n if(nops( punt)<>0) then\n plot_arcering:=plots[listplot]([seq(punt[n],n=1. .k-1)],\n style=point,color=green);\n end if; \n plot_nietarcering:=plots[listplot]([seq(nietpunt[n],n=1..l-1)], \n style=point,color=red);\n if(nops(punt)=0 ) then\n plots[display](\{plot_gebied,plot_nietarcering\});\n \+ else \n plots[display](\{plot_gebied,plot_arcering,plot_nietarcer ing\});\n end if;\n end if;\nend proc:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 875 "stabintervalPC:=proc(PCpaar) \nlocal vt, grens, st ep, j, ok, opl, tol, L_,mu_;\n mu_:=PCpaar['mu'];\n L_:=PCpaar['L']; \n if ((L_<>\"Lin\") and (L_<>\"n\") and (L_<>\"Lout\")) then\n ER ROR(\"De laatste parameter in PCpaar kan enkel \\\"n\\\", \\\"Lin\\\" \+ of \\\"Lout\\\" zijn.\");\n end if; \n if ((mu_=0) and (L_<>\"n\")) \+ then\n ERROR(\"Men kan enkel lokale extrapolatie toepassen als men \+ minstens \351\351nmaal de corrector toepa st\");\n end if;\n vt:=stabvtPC(PCpaar);\n tol:=1/10.^12;\n step:= 0.2;\n j:=0;\n ok:=true;\n grens:=0;\n while step>tol do\n whil e ok do\n j:=j+1; \n opl:=solve(subs(h=-(grens+j*step),vt)=0 ,r);\n if sort(map(abs,[opl]))[-1]>1 then \n ok:=false; \n end if;\n end do;\n grens:=grens+(j-1)*step;\n j:=0;\n \+ ok:=true;\n step:=step/10.;\n end do; \n printf(`stabiliteitsi nterval :`);\n return([-grens,0]);\nend proc:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 1748 "helpPC:=proc()\nprint(`Dezelfde afspraken i. v.m. de input van methodes blijft gelden (typ \\`helpLMM( )\\` voor me er uitleg).\\n`);\nprint(`Besproken procedures: \\`ordePC\\`, \\`stabv tPC\\`, \\`stabgebiedPC\\` ,\\`randkrommePC\\` en \\`stabintervalPC\\` .\\n`);\nprint(`De procedure \\`ordePC\\` berekent de orde van het PC- paar. Als input geef je een lijst met 5 argumenten, nl. \\`P\\`, \\`C \\`, \\`mu\\`, \\`t\\` en \\`L\\`. Men kan dit op de volgende manier d oen:\\n`);\nprintf(` PCpaar['P']:=AdamsBashforth(p=2);\\n`);\nprin tf(` PCpaar['C']:=AdamsMoulton(k=4);\\n`);\nprintf(` PCpaar['m u']:=3;\\n`);\nprintf(` PCpaar['t']:=0;\\n`);\nprintf(` PCpaar ['L']:=\"Lin\";\\n`);\nprint(`Let op: de namen \\`P\\`, \\`C\\`, \\`mu \\`, \\`t\\` en \\`L\\` mogen niet veranderd worden en de enkele quote s errond zijn ook belangrijk! Om de procedure te gebruiken kan je dan \+ ingeven:\\n`); printf(` ordePC(PCpaar);`);\nprint(`Voor deze L zi jn er drie mogelijkheden:\\n`);\nprint(`\\\"n\\\": geen lokale extrapo latie, `);\nprint(`\\\"Lin\\ \": in elke loop gebeurt lokale extrapolatie, `);\nprint (`\\\"Lout\\\": enkel lokale extrapolatie op het einde, na de lus.`); \nprint(`Voor de procedures \\`stabvtPC\\`, \\`stabgebiedPC\\`, \\`ran dkrommePC\\` en \\`stabintervalPC\\` geldt dezelfde afspraak. Als inpu t bij \\`stabgebiedPC\\` is er nog een optionele parameter \\`GRENSWAA RDE\\` analoog met \\`stabgebied\\` (zie helpLMM( ) ). \\`stabgebiedPC \\` is een uitgebreide procedure en kan dus enige tijd in beslag nemen .\\n`);\nprintf(` ordePC:=proc(PCpaar)\\n`);\nprintf(` stabvtP C:=proc(PCpaar)\\n`);\nprintf(` stabgebiedPC:=proc(PCpaar,GRENSWAA RDE::realcons)\\n`);\nprintf(` randkrommePC:=proc(PCpaar)\\n`);\np rintf(` stabintervalPC:=proc(PCpaar)\\n`);\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2676 "take_stepPC:=proc(n::posint,h_,pr obleem,PCpaar)\n local opl,f,W,y,k,i,j,f_j,take_stepLMM,f_t,k_,y_j,l, A,x1,y1,P_,C_,mu_,t_,L_;\n A:=probleem['matrix'];\n x1:=probleem['x0 '];\n y1:=probleem['y0'];\n P_:=PCpaar['P'];\n C_:=PCpaar['C'];\n \+ mu_:=PCpaar['mu'];\n t_:=PCpaar['t'];\n L_:=PCpaar['L'];\n if ((L_< >\"n\") and (L_<>\"Lin\") and (L_<>\"Lout\")) then\n ERROR(\"De laa tste parameter in PCpaar mag enkel \\\"n\\\", \\\"Lin\\\" of \\\"Lout \\\" zijn.\");\n end if;\n if ((mu_=0) and (L_<>\"n\")) then\n ER ROR(\"Men kan enkel lokale extrapolatie toepassen als men minstens \+ \351\351nmaal de corrector toepast.\");\n end if;\n k _:=linalg[coldim](P_)-1;\n k:=linalg[coldim](C_)-1;\n if (P_[2,k_+1] <>0) then\n ERROR(\"De predictor moet een expliciete methode zijn\" );\n elif (C_[2,k+1]=0) then\n ERROR(\"De corrector moet een impli ciete methode zijn\");\n end if;\n f:=(y)->linalg[multiply](probleem ['matrix'],y);\n W:=foutconstante(C_)/(foutconstante(P_)-foutconstant e(C_));\n\n y:=genereer_beginwaarden(P_,probleem,h_);\n for i from 1 to k_ do\n f_t[i]:=f(y[i]);\n end do;\n if (k>k_) then\n opl: =exact(probleem,x);\n for j from k_+1 to k do\n y:=[seq(y[i],i =1..j-1),eval(subs(\{x=x1+(j-1)*h_\},opl))];\n f_t[j]:=f(y[j]);\n end do;\n end if;\n\n take_stepLMM:=proc(n::integer,h_,A,x1,y1,m ethode::matrix,y,f_t)\n local expr,i,k_methode;\n k_methode:=lin alg[coldim](methode)-1;\n expr:=array(sparse,1..linalg[coldim](A)); \n for i from 1 to k_methode do\n expr:=linalg[matadd](expr,li nalg[matadd](methode[1,i]*y[n-k_methode+i-1],\n h_*methode[ 2,i]*f_t[n-k_methode+i-1],-1,1));\n end do;\n return([seq(y[i],i =1..n-1),evalm(expr)]);\n end proc:\n\n if (mu=0) then\n for j fr om k_+1 to n do\n y:=take_step(j,P_,probleem,h_,y);\n end do; \n return(y);\n end if;\n for j from max(k_,k)+1 to n do\n y_j [0]:=take_stepLMM(j,h_,A,x1,y1,P_,y,f_t)[j];\n for i from 1 to mu_ \+ do\n f_j[i-1]:=f(y_j[i-1]);\n y_j[i]:=array(sparse,1..linalg [coldim](probleem['matrix']));\n for l from 1 to k do\n y_ j[i]:=linalg[matadd](y_j[i],linalg[matadd](C_[1,l]*y[j-k+l-1],h_*C_[2, l]*f_t[j-k+l-1],-1,1));\n end do;\n y_j[i]:=linalg[matadd](y _j[i],C_[2,k+1]*h_*f_j[i-1]);\n\n if (L_ = \"Lin\") then\n \+ y_j[i]:=linalg[matadd]((1+W)*y_j[i],W*y_j[0],1,-1);\n end if;\n \+ end do;\n if (L_=\"Lout\") then\n y_j[mu_]:=linalg[matadd]( (1+W)*y_j[mu_],W*y_j[0],1,-1); \n end if;\n y:=[seq(y[i],i=1. .j-1),eval(y_j[mu_])];\n if (t_=1) then\n f_t[j]:=f_j[mu_-1]; \n elif (t_=0) then\n f_t[j]:=f(y[j]);\n else\n ERROR( `t (vierde parameter in PCpaar) kan enkel 0 of 1 zijn`);\n end if; \n end do;\n RETURN(eval(y));\nend proc:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 398 "LAFPC:=proc(probleem,PCpaar,m::posint)\nlocal y,k, l,opl,x1,P_,C_;\nx1:=probleem['x0'];\nP_:=PCpaar['P'];\nC_:=PCpaar['C' ];\nif (nargs=3) then\n l:=m;\nelse\n l:=max(orde(P_),orde(C_))+3;\n end if;\nk:=max(linalg[coldim](P_),linalg[coldim](C_))-1;\ny:=take_ste pPC(k+1,'h',probleem,PCpaar);\nopl:=exact(probleem,x);\nreturn(map(x-> series(x,h=0,l),\n evalm(array([op(subs(x=x1+(k)*'h',opl))])-y[k+1])) );\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 423 "numPC:=pr oc(np::posint,h::realcons,probleem,PCpaar)\nlocal y,n,k,i,plot,A,x1,y1 ,P_,C_,mu_,t_,L_;\nA:=probleem['matrix'];\ny:=take_stepPC(np,h,problee m,PCpaar);\nfor i from 1 to linalg[coldim](A) do\nplot[i]:=plots[listp lot]([seq([j*h,y[j][i]],j=1..np)],style=point,color=COLOR(RGB, rand()/ 10^12, rand()/10^12, rand()/10^12)); \nen d do;\nplots[display](\{seq(plot[i],i=1..linalg[coldim](A))\});\nend p roc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 508 "GAFPC:=proc(np::po sint,h::realcons,probleem,PCpaar)\n local y, opl,i,plot,A,x1;\n A:=p robleem['matrix'];\n x1:=probleem['x0'];\n y:=take_stepPC(np,h,probl eem,PCpaar);\n opl:=exact(probleem,x);\n for i from 1 to linalg[cold im](A) do\n plot[i]:=plots[listplot]([seq([j*h,subs(x=x1+j*h,opl[i] )-y[j+1][i]],j=1..np-1)],\n style=point,color=COLOR(RGB, r and()/10^12, rand()/10^12, rand()/10^12)); \+ end do;\n plots[display](\{seq(plot[i],i=1..linalg[coldim](A))\}) ;\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 3019 "opvolging :=proc(PCpaar,m::integer)\nlocal l,W,A,ordeP,ordeC,k,x1,y1,f,j,y,opl,t ake_stepLMM,i,f_j,y_j,q,f_t,s,k_,P_,C_,mu_,t_,L_,PC,probleem;\nP_:=PCp aar['P'];\nC_:=PCpaar['C'];\nmu_:=PCpaar['mu'];\nt_:=PCpaar['t'];\nL_: =PCpaar['L'];\nif ((t_<>0) and (t_<>1)) then\n ERROR(\"t (voorlaatste parameter in PCpaar) kan enkel 0 of 1 zijn.\");\nend if;\nif (nargs=2 ) then\n l:=m;\nelse\n l:=max(orde(P_),orde(C_))+3;\nend if;\nproble em['matrix']:=linalg[matrix]([[1]]);\nprobleem['y0']:=[1];\nprobleem[' x0']:=0;\nf:=(y)->linalg[multiply](A,y);\nordeP:=orde(P_);\nordeC:=ord e(C_);\nprintf(`predictor: (orde %d)`,ordeP);\nprint(P_);\nprintf(`cor rector: (orde %d)`,ordeC);\nprint(C_);\nif (mu_<>0) then\n if (L_=\"n \") then\n if (t_=1) then\n printf(`uitgevoerde mode: P(EC)^%d \\n\\n`,mu_);\n else\n printf(`uitgevoerde mode: P(EC)^%dE\\n \\n`,mu_);\n end if;\n elif (L_=\"Lin\") then\n if (t_=1) then \n printf(`uitgevoerde mode: P(ECL)^%d\\n\\n`,mu_);\n else\n \+ printf(`uitgevoerde mode: P(ECL)^%dE\\n\\n`,mu_);\n end if;\n \+ elif (L_=\"Lout\") then\n if (t_=1) then\n printf(`uitgevoerde mode: P(EC)^%dL\\n\\n`,mu_);\n else\n printf(`uitgevoerde mod e: P(EC)^%dLE\\n\\n`,mu_);\n end if;\n else\n ERROR(\"De laatst e parameter kan enkel \\\"n\\\", \\\"Lin\\\" of \\\"Lout\\\" zijn.\"); \n end if;\nelse\n if ((L_=\"n\") or (L_=\"Lin\")) then\n printf( `uitgevoerde mode: P\\n\\n`);\n elif (L=\"Lout\") then\n ERROR(`De mode \\`PL\\` bestaat niet`);\n else\n ERROR(\"De laatste paramet er kan enkel \\\"n\\\", \\\"Lin\\\" of \\\"Lout\\\" zijn.\");\n end i f;\nend if;\nprintf(`Fout toegepast op het probleem y'=y:\\n`);\nk:=ma x(linalg[coldim](P_),linalg[coldim](C_));\nopl:=x->exp(x);\nprintf(`na de predictor:`);\nprint((LAF(P_,probleem,l)[1]));\nif (L_<>\"Lin\") t hen\n if (mu_<>0) then\n for j from 1 to mu_ do\n printf(`ite ratie %d:\\n`,j);\n printf(` fout na C:\\n`);\n PC['P']:=P_ ; PC['C']:=C_; PC['mu']:=j; PC['t']:=t_; PC['L']:=\"n\";\n print( LAFPC(probleem,PC,l)[1]);\n end do;\n if (L_=\"Lout\") then\n \+ printf(`Lokale extrapolatie:\\n`);\n printf(` fout na L:\\n`) ;\n PC['P']:=P_; PC['C']:=C_; PC['mu']:=mu_; PC['t']:=t_; PC['L'] :=\"Lout\";\n print(LAFPC(probleem,PC,l)[1]);\n end if;\n end if;\nelse\n if (mu_<>0) then\n PC['P']:=P_; PC['C']:=C_; PC['mu'] :=1; PC['t']:=1; PC['L']:=\"n\";\n y:=take_stepPC(k,'h',probleem,PC );\n for s from 1 to mu_ do\n printf(`iteratie %d:\\n`,s);\n \+ printf(` fout na C:\\n`);\n print(simplify(map(x->series(x,h =0,l),\n evalm(linalg[matadd]([opl(probleem['x0']+(k-1 )*'h')],-y[k]))))[1]);\n printf(` fout na L:\\n`);\n PC['P' ]:=P_; PC['C']:=C_; PC['mu']:=s; PC['t']:=1; PC['L']:=\"Lin\";\n \+ print(LAFPC(probleem,PC,l)[1]);\n y:=take_stepPC(k,'h',probleem,P C);\n k_:=linalg[coldim](C_)-1;\n y_j:='h'*C_[2,k_+1]*[y[k][ 1]];\n for q from 1 to k_ do\n y_j:=linalg[matadd](y_j,lin alg[matadd](C_[1,q]*y[k-k_+q-1],'h'*C_[2,q]*y[k-k_+q-1],-1,1));\n \+ end do;\n y:=[seq(y[i],i=1..k-1),eval(y_j)];\n end do;\n end if;\nend if;\nreturn(NULL);\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1499 "helpOP_PROBLEEM_PC:=proc()\nprint(`Dezelfde afsprak en i.v.m. de input van methodes blijft gelden (typ \\`helpLMM( )\\` vo or meer uitleg), voor de input van problemen (typ \\`helpOP_PROBLEEM_L MM( )\\` voor meer uitleg) en voor de input van PC-paren (typ \\`helpP C( )\\`voor meer uitleg).\\n`);\nprint(`Besproken procedures: \\`take_ stepPC\\`, \\`LAFPC\\`, \\`numPC\\`, \\`GAFPC\\` en \\`opvolging\\`.\\ n`);\nprint(`De procedure \\`take_stepPC\\` berekent alle numerieke wa arden tot een gegeven punt (analoog met \\`take_step2\\` uit \\`LMM.m \\`). Als input geef je het aantal te berekenen punten, de \\`h\\`-waa rde, het probleem en het PC-paar. Aan \\`LAFPC\\` (lokale afknottingsf out) geef je het probleem, het PC-paar en (optioneel) een getal tot wa ar je de LAF wil ontwikkelen. \\`numPC\\` geeft een plot van de numeri eke oplossing en \\`GAFPC\\` van de globale afknottingsfout. Telkens g eef je het aantal punten, de \\`h\\`-waarde, het probleem en het PC-pa ar. Opvolging is een procedure die per stap dat er uitgevoerd wordt, d e lokale afknottingsfout op het probleem y\\' = y berekent. Deze proce dure heeft enkel het PCpaar nodig, en optioneel een getal tot waar de \+ ontwikkeling van de LAF moet gebeuren.`);\nprintf(` take_stepPC:=p roc(n::posint,h_,probleem,PCpaar)\\n`);\nprintf(` LAFPC:=proc(prob leem,PCpaar,m::posint)\\n`);\nprintf(` numPC:=proc(np::posint,h::r ealcons,probleem,PCpaar)\\n`);\nprintf(` GAFPC:=proc(np::posint,h: :realcons,probleem,PCpaar)\\n`);\nprintf(` opvolging:=proc(PCpaar, m::integer)\\n`);\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1300 "schur:=proc(poly::polynom,x::symbol)\nlocal constructmatrix,re curs,dim,polynoom;\nglobal lst;\n\nconstructmatrix:=proc(poly,x,a)\n \+ local subspoly,z,k,r0,j,i,l;\n subspoly:=subs(x=(1+z)/(1-z),poly);\n \+ k:=degree(poly,x);\n subspoly:=expand(simplify((1-z)^k*subspoly)); \+ \n r0:=coeff(subspoly,z,k);\n if type(r0,numeric) then\n if r0<0 \+ then \n subspoly:=-subspoly;\n end if;\n end if;\n a:=array( 1..k,1..k);\n for i from 1 to k do \n for j from 1 to k do \n \+ l:=i+(j-i)*2;\n if (l*(k-l)<0) then \n a[i,j]:=0\n e lse \n a[i,j]:=coeff(subspoly,z,k-l);\n end if\n end do ;\n end do;;\n RETURN();\nend proc:\n\nrecurs:= proc(mat::matrix,k:: posint)\n local l,hulpmat;\n global lst;\n if k=1 then \n lst:=l st union \{factor(mat[1,1])\}\n else\n for l from 1 to k do\n \+ lst:=lst union \{factor(linalg[det](linalg[minor](mat,l,l)))\};\n \+ hulpmat:=linalg[delcols](linalg[delrows](mat,l..l),l..l);\n rec urs(hulpmat,k-1);\n end do;\n end if;\n RETURN();\nend proc:\n\np olynoom:=expand(poly);\nif degree(polynoom,x)>1 then\n constructmatri x(poly,x,mat);\n print(mat);\n lst:=\{factor(mat[2,1])\};\n dim:=no ps(convert(mat,listlist));\n recurs(mat,dim);\n lst:=map(u->factor(u >=0),lst);\n RETURN(lst);\nelse \n RETURN(\{-1 <=solve(polynoom,x),s olve(polynoom,x)<=1\});\nend if:\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 169 "uitschrijven:=proc(methode::matrix)\nlocal k;\nk: =linalg[coldim](methode);\nreturn(sum(methode[1,k+1-i]*y[n-i+2],i=1..k )=h*sum(methode[2,k+1-i]*f[n-i+2],i=1..k));\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 515 "save uitschrijven,G, gam, beta, Ad amsBashforth, AdamsMoulton, Nystrom, MilneSimpson, Maximaal, uitschrij ven, orde, randkromme, plotten, stabgebied, stabvt, controleer, PLAF, \+ foutconstante, BDF, stabinterval, nulstabiel, helpLMM, exact, genereer _beginwaarden, take_step, take_step2, num, LAF, GAF, helpOP_PROBLEEM_L MM ,rho, sigma, zetstapgelijk, simple, stabvtPC, ordePC, stabgebiedPC, randkrommePC, stabintervalPC, helpPC, take_stepPC, LAFPC, numPC, GAFP C, opvolging, helpOP_PROBLEEM_PC,schur, \"C:\\\\temp\\\\LiMMap.m\";" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "45 0 0" 18 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }