|
1 | 1 | #randomgrowth2.jl
|
2 | 2 | using PyPlot
|
3 |
| -using ODE |
| 3 | +using OrdinaryDiffEq |
4 | 4 |
|
5 | 5 | function randomgrowth2()
|
6 | 6 | num_trials = 2000
|
@@ -32,16 +32,22 @@ function randomgrowth2()
|
32 | 32 | t0 = 4
|
33 | 33 | tn = -6
|
34 | 34 | dx = 0.005
|
35 |
| - deq = (t, y) -> [y[2]; t*y[1]+2*y[1]^3; y[4]; y[1]^2] |
36 |
| - y0 = [airy(t0); airy(1,t0); 0; airy(t0)^2] # boundary conditions |
37 |
| - t, y = ode23(deq, y0, t0:-dx:tn) # solve |
38 |
| - F2 = Float64[exp(-y[i][3]) for i = 1:length(y)] # the distribution |
39 |
| - f2 = gradient(F2, t) # the density |
| 35 | + deq = function (t, y, dy) |
| 36 | + dy[1] = y[2] |
| 37 | + dy[2] = t*y[1]+2y[1]^3 |
| 38 | + dy[3] = y[4] |
| 39 | + dy[4] = y[1]^2 |
| 40 | + end |
| 41 | + y0 = big.([airy(t0); airy(1,t0); 0; airy(t0)^2]) # boundary conditions |
| 42 | + prob = ODEProblem(deq,y0,(t0,tn)) |
| 43 | + sol = solve(prob,Vern8(), saveat=-dx, abstol=1e-12, reltol=1e-12) # solve |
| 44 | + F2 = Float64[exp(-sol[i][3]) for i = 1:length(y)] # the distribution |
| 45 | + f2 = gradient(F2, t) # the density |
40 | 46 |
|
41 | 47 | # add(p, Curve(t, f2, "color", "red", "linewidth", 3))
|
42 | 48 | # Winston.display(p)
|
43 | 49 | subplot(1,length(Ns),jj)
|
44 |
| - plot(t, f2, "r", linewidth = 3) |
| 50 | + plot(sol.t, f2, "r", linewidth = 3) |
45 | 51 | ylim(0, 0.6)
|
46 | 52 | println(mean(C))
|
47 | 53 | end
|
|
0 commit comments