> #Schiessverfahren
> restart;
> with(linalg): with(plots):
> Schiessverfahren:=proc(a,b,f,fy,fyy,c1,c2,s,N,anz)
> # Löse y''(x)=f(x,y,y') auf Intervall (a,b)
> # f=f(x,y,yy): R^3-> R (yy ist die erste Ableitung y')
> # fy=fy(x,y,yy) Ableitung von f nach y
> # fyy=fyy(x,y,yy) Ableitung von f nach yy=y'
> # Randbedingungen y(a)=c1,y(b)=c2
> # s erste Schätzung des Anfangswertes s=y'(a)
> # Schrittweite ist (b-a)/N
> # anz Anzahl der Schiessversuche
> local RungeKuttasystem,k,j,w0,f1,w,w1,z,h;
> global p,s0,s1;
> # option trace;
> RungeKuttasystem:=proc(y0,t0,f,T,N)
> # Schrittweite (T-t0)/N
> # option trace;
> local d,k,y,t,k1,k2,k3,k4,j,h;
> global p;
> d:=vectdim(y0); p:=vector(N+1); p[1]:=y0[1];
> y:=vector(d); 
> k1:=vector(d); k2:=vector(d); k3:=vector(d); k4:=vector(d); 
> for j from 1 to d do y[j]:=y0[j] od;
> t:=t0; h:=(T-t0)/N;
> for j from 0 to N-1 do
>  k1:=f(t,y); 
>  k2:=f(t+h/2,matadd(y,k1,1,h/2)); 
>  k3:=f(t+h/2,matadd(y,k2,1,h/2));
>  t:=t+h; 
>  k4:=f(t,matadd(y,k3,1,h)); 
>  for k from 1 to d do y[k]:=y[k]+h*(k1[k]+2*k2[k]+2*k3[k]+k4[k])/6 od;
>  p[j+2]:=y[1]
> od;
> evalm(y)
> end:
> h:=(b-a)/N;
> p:=vector(N+1); z:=vector(N+1);
> # Erste Schätzung des Anfangswertes s=y'(a)
> s1:=s; print(1,`-te Schätzung von s=y'(a):`,s1);
> for k from 1 to anz do
>  # Aufstellen des AWP-Systems
>  w0:=vector(4);
>  w0[1]:=c1; w0[2]:=s1; w0[3]:=0; w0[4]:=1;
>  f1:=vector(4); w:=vector(4);
>  f1:=proc(x,w) local g; 
>  g:=vector(4); 
>  g[1]:=w[2]; g[2]:=f(x,w[1],w[2]);
>  g[3]:=w[4]; g[4]:=fy(x,w[1],w[2])*w[3]+fyy(x,w[1],w[2])*w[4];
>  evalm(g) end;
>  # Berechnung des AWP-Systems
>  w1:=vector(4);
>  w1:=RungeKuttasystem(w0,a,f1,b,N);
>  s0:=s1;
>  #neue Schätzung des Anfangswertes s=y'(a)
>  s1:=s1-(w1[1]-c2)/w1[3]; print(k+1,`-te Schätzung von s=y'(a):`,s1);
> od;
> #graphische Ausgabe
> #print(p);
> print(`erreichter Randwert y(b) für`,anz,`te Schätzung`, s0,`:`,p[N+1]);
> for j from 1 to N+1 do z[j]:=a+(j-1)*h od;
> listplot([seq([z[j],p[j]],j=1..N+1)],thickness=0);
> end:
Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

> #Beispiel 1
> # -y''(x)+y(x)*y'(x)=-1 in (0,1), y(0)=y(1)=0
> a:=0.0: b:=1.0:
> f:=proc(x,y,yy) local g; g:=-1+y*yy; RETURN(g) end:
> fy:=proc(x,y,yy) local g; g:=-yy; RETURN(g) end:
> fyy:=proc(x,y,yy) local g; g:=-y; RETURN(g) end:
> c1:=0.0: c2:=0.0: s:=0.3: N:=20: 
> Schiessverfahren(a,b,f,fy,fyy,c1,c2,s,N,1);
> Schiessverfahren(a,b,f,fy,fyy,c1,c2,s,N,3);
> 

                  1, -te Schätzung von s=y'(a):, .3


              2, -te Schätzung von s=y'(a):, .4927151105


  erreichter Randwert y(b) für, 1, te Schätzung, .3, :, -.1975369276



                  1, -te Schätzung von s=y'(a):, .3


              2, -te Schätzung von s=y'(a):, .4927151105


              3, -te Schätzung von s=y'(a):, .4961003360


              4, -te Schätzung von s=y'(a):, .4958108665


  erreichter Randwert y(b) für, 3, te Schätzung, .4961003360, :,

        .00027772131


> 
> 
>