|
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