Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- (*Mathematica*)(*start*)Clear[A, B, a, b, x];
- a = {1, 1, 1, 1, -3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0, 0};
- b = (Range[Length[a]] - 1)!;
- a = a*b;
- nn = Length[a];
- d = Max[Flatten[Position[Sign[Abs[a]], 1]]]
- A = Table[
- Table[If[n >= k, Binomial[n - 1, k - 1]*a[[n - k + 1]], 0], {k, 1,
- Length[a]}], {n, 1, Length[a]}];
- b = Inverse[A][[All, 1]];
- x = N[Table[(n - 1)*b[[n - 1]]/b[[n]], {n, nn - 8, nn - 1}], 30]
- perturbation = -1;
- x - perturbation
- Sum[a[[k + 1]]/k!*x^k, {k, 0, d}]
- Clear[x];
- polynomial = Sum[a[[k + 1]]/k!*x^k, {k, 0, d}]
- a1 = CoefficientList[
- Expand[Sum[a[[k + 1]]/k!*(x + perturbation)^k, {k, 0, d}]],
- x]*(Range[d] - 1)!
- Clear[A, B, a, b, x];
- a = Flatten[{a1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0}];
- nn = Length[a];
- d = Max[Flatten[Position[Sign[Abs[a]], 1]]]
- A = Table[
- Table[If[n >= k, Binomial[n - 1, k - 1]*a[[n - k + 1]], 0], {k, 1,
- Length[a]}], {n, 1, Length[a]}];
- b = Inverse[A][[All, 1]];
- x = N[Table[(n - 1)*b[[n - 1]]/b[[n]], {n, nn - 8, nn - 1}], 30]
- "perturbed x polynomial"
- Sum[a[[k + 1]]/k!*x^k, {k, 0, d}]
- x = x + perturbation
- "unperturbed x polynomial"
- polynomial
- Clear[x];
- Sum[a[[k + 1]]/k!*x^k, {k, 0, d}]
- (*end*)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement