CA II - Uebungsblatt 07

Aufgabe 1

(a)

f=(y[x]-r)^2+x^2-r^2

rule=(y[x]^s_/;s≥2)->Expand[ -(f-y[x]^2)y[x]^(s-2)]

k=1;
r[k]=Solve[D[f,{x,k}]==0,D[y[x],{x,k}]][[1,1]]

Ansatz=D[y[x],{x,k}]+Sum[a[j]D[y[x],{x,j}],{j,0,k-1}]

Poly=Numerator[Together[Ansatz/.Table[r[j],{j,k}]]]

Poly=Poly//.rule

Solve[Map[#==0&,CoefficientList[Poly,y[x]]],Table[a[j],{j,0,k-1}]]

Naechster Versuch mit k=2.

k=2;
r[k]=Solve[D[f,{x,k}]==0,D[y[x],{x,k}]][[1,1]]

Ansatz=D[y[x],{x,k}]+Sum[a[j]D[y[x],{x,j}],{j,0,k-1}]

Poly=Numerator[Together[Ansatz//.Table[r[j],{j,k}]]]

Poly=Poly//.rule

Solve[Map[#==0&,CoefficientList[Poly,y[x]]],Table[a[j],{j,0,k-1}]][[1]]

DE=Collect[Numerator[Together[Ansatz/.%]],D[y[x],{x,_}],Factor]==0

DEtoRE[DE,y[x],B[m]]

Die Rekursion ist vom hypergeometrischen Typ mit Symmetriezahl 2.

(b)

f=y[x]^3+x y[x]^2+x^2

rule=(y[x]^s_/;s≥3)->Expand[ -(f-y[x]^3)y[x]^(s-3)]

Fuer alle k<3 ist die Suche nach einer Differentialgleichung nicht erfolgreich!

k=1;
r[k]=Solve[D[f,{x,k}]==0,D[y[x],{x,k}]][[1,1]]
k=2;
r[k]=Solve[D[f,{x,k}]==0,D[y[x],{x,k}]][[1,1]]
k=3;
r[k]=Solve[D[f,{x,k}]==0,D[y[x],{x,k}]][[1,1]]

Ansatz=D[y[x],{x,k}]+Sum[a[j]D[y[x],{x,j}],{j,0,k-1}]

Poly=Numerator[Together[Ansatz//.Table[r[j],{j,k}]]]

Poly=Poly//.rule;

Solve[Map[#==0&,CoefficientList[Poly,y[x]]],Table[a[j],{j,0,k-1}]][[1]]

DE=Collect[Numerator[Together[Ansatz/.%]],D[y[x],{x,_}],Factor]==0

DEtoRE[DE,y[x],B[m]]

Diese Rekursion ist sogar von erster Ordnung, also insbesondere vom hypergeometrischen Typ.

(c)

Exit[]

Needs["SpecialFunctions`"]

AlgebraicDE[(y[x]-r)^2+x^2-r^2,y[x]]

AlgebraicDE[y[x]^3+x y[x]^2+x^2,y[x]]

Aufgabe 2

(a)

Clear[f];
f[x_]:=x^3-7

Df[x_]:=D[f[y],y]/.y->x

L=FixedPointList[#-f[#]/Df[#]&,N[2,20]]

Table[7^(1/3),{6}]-L

Die quadratische Konvergenz erkennt man daran, dass sich der Fehler bei jedem Iterationsschritt um das Quadrat verkleinert.

(b)

F=x^2 y+y^2 x-y-1

(i)

sol=Solve[F==0,y]

Fexplicit1=y/.sol[[1]]

Plot[Fexplicit1,{x,0,4}]

Fexplicit2=y/.sol[[2]]

Plot[Fexplicit2,{x,0,4}]

Als Loesung kommt wegen y[0]=-1 nur die erste explizite Loesung in Frage.

(ii)

Normal[Series[Fexplicit1,{x,0,15}]]

benoetigte Funktion

MAXDEORDER=5

Clear[HolonomicDE];
HolonomicDE[f_,F_[x_],print_:False]:=Module[{DF,eqcl,rat,A,approach,sol,eq,i,j,l,m},
DF=Table[D[f,{x,i}],{i,0,MAXDEORDER-1}];
DF=Map[If[Head[#]===Plus,Apply[List,#],{#}]&,DF];(* DF consists of lists of derivatives (sums are lists of summands) *)
If[Head[f]==Plus,eqcl=Table[{A[1] f[[i]]},{i,Length[f]}],eqcl={A[1]f}];
For[i=2,i≤MAXDEORDER,i++,
For[l=1,l≤Length[DF[[i]]],l++, (* for every summand of the (i-1)-th derivative do the following *)
m=Length[eqcl];
For[j=1,j≤m,j++,
rat=Together[DF[[i,l]]/eqcl[[j,1]]]; (* building quotient of l-th summand and first element of j-th eq. class*)
If [PolynomialQ[Denominator[rat],x] && PolynomialQ[Numerator[rat],x], (* test, if rat is rational or not  *)
AppendTo[eqcl[[j]],A[i]DF[[i,l]]]; (* add to equivalence class *)
Break[],
If[j==Length[eqcl],AppendTo[eqcl,{A[i]DF[[i,l]]}];  (* build new equivalence class *)
];
];
];
];
var=Table[A[k],{k,i}];
approach=Map[Apply[Plus,#]&,eqcl]; (* sum over all elements in equivalence classes *)
approach=Map[#==0&,approach];
If[print,Print[{eqcl,approach}]];
Off[Solve::"svars"];
sol=Solve[approach,var];
If[sol=!={{}} && Union[var/.sol[[1]]]=!={0},
eq=Sum[A[k+1]D[F[x],{x,k}],{k,0,i-1}]/.sol[[1]];
eq=Numerator[Together[eq]]/.A[_]->1;
Return[Collect[eq,D[F[x],{x,_}],Factor]==0]];
];
]

(iii)

HolonomicDE[(1-x^2)/(2x)-Sqrt[x^4-2x^2+4x+1]/(2x),G[x],True]

(iv)

F=F/.y->y[x]

rule=(y[x]^s_/;s≥2)->Expand[ -(F-x*y[x]^2)/x y[x]^(s-2)]

k=1;
r[k]=Solve[D[F,{x,k}]==0,D[y[x],{x,k}]][[1,1]]
k=2;
r[k]=Solve[D[F,{x,k}]==0,D[y[x],{x,k}]][[1,1]]

Ansatz=D[y[x],{x,k}]+Sum[a[j]D[y[x],{x,j}],{j,0,k-1}]

Poly=Numerator[Together[Ansatz//.Table[r[j],{j,k}]]]

Poly=Poly//.rule;

Solve[Map[#==0&,CoefficientList[Poly,y[x]]],Table[a[j],{j,0,k-1}]][[1]]

DE=Collect[Numerator[Together[Ansatz/.%]],D[y[x],{x,_}],Factor]==0

benoetigte Funktionen

Clear[ImplicitTaylor];
ImplicitTaylor[F_,y_[x_],y0_,n_]:=Module[{ystrich,k},
ystrich=-D[F,x]/D[F,y];
ableitungen=NestList[Together[D[#,x]+D[#,y]*ystrich]&,ystrich,n-1];
y0+Sum[Limit[Limit[ableitungen[[k]],y->y0],x->0]/k!*x^k,{k,1,n}]
];

Clear[FastImplicitTaylor];
FastImplicitTaylor[f_,y_[x_],y0_,n_]:=Module[{F,z,approx},
F[z_]:=Normal[(f/.y->z)+O[x]^(n+1)];
approx=y0;
Do[approx=Normal[#-F[#]/F'[#]&[approx]+O[x]^(2^k)],
{k,1,Log[2,n+1]}];
approx+O[x]^(n+1)
]/;IntegerQ[Log[2,n+1]]

(v)

Timing[ImplicitTaylor[F,y[x],-1,25]][[1]]

Timing[FastImplicitTaylor[F,y[x],-1,200]][[1]]