solvers/NewtonImplicitSolver.py

This example looks at a use of rootfinders

View output (PDF) | source (python)

See also
casadi::rootfinder();
casadi::rootfinder
CASADI_EXPORT Function rootfinder(const std::string &name, const std::string &solver, const SXDict &rfp, const Dict &opts=Dict())
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 from casadi import *
28 from numpy import *
29 from pylab import *
30 
31 
32 #!
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 eps = SX.sym("eps")
44 mu = SX.sym("mu")
45 alpha = SX.sym("alpha")
46 k = SX.sym("k")
47 sigma = SX.sym("sigma")
48 params = [eps,mu,alpha,k,sigma]
49 
50 
51 
52 a = SX.sym("a")
53 gamma = SX.sym("gamma")
54 
55 
56 
57 res0 = mu*a+1.0/2*k*a*sin(gamma)
58 res1 = -sigma * a + 3.0/4*alpha*a**3+k*a*cos(gamma)
59 
60 
61 
62 sigma_ = 0.1
63 alpha_ = 0.1
64 k_ = 0.2
65 params_ = [0.1,0.1,alpha_,k_,sigma_]
66 
67 
68 
69 f=Function("f", [vertcat(a,gamma),vertcat(*params)],[vertcat(res0,res1)])
70 opts = {}
71 opts["abstol"] = 1e-14
72 opts["linear_solver"] = "csparse"
73 s=rootfinder("s", "newton", f, opts)
74 
75 
76 
77 x_ = s([1,-1], params_)
78 
79 print("Solution = ", x_)
80 
81 
82 
83 x = [sqrt(4.0/3*sigma_/alpha_),-0.5*pi]
84 print("Reference solution = ", x)
85 
86 
87 
88 residual = f(x_, params_)
89 print("residual = ", residual)
90 
91 for i in range(1):
92  assert(abs(x_[i]-x[i])<1e-6)
93 
94 
95 
96 print(s.stats())