1 Julia

1.2 Connect Julia

library(JuliaConnectoR)
library(JuliaCall)
              _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.7.0-beta3.0 (2021-07-07)
 _/ |\__'_|_|_|\__'_|  |  release-1.7/e76c9dad42 (fork: 68 commits, 62 days)
|__/                   |

julia>

2 Demontration: Run Julia in Rstudio

2.1 Basic Plot

using Plots
plot(Plots.fakedata(50, 5), w = 3)

2.2 Marginal Histogram

using RDatasets, StatsPlots
iris = dataset("datasets", "iris")
150×5 DataFrame
 Row │ SepalLength  SepalWidth  PetalLength  PetalWidth  Species
     │ Float64      Float64     Float64      Float64     Cat…
─────┼─────────────────────────────────────────────────────────────
   1 │         5.1         3.5          1.4         0.2  setosa
   2 │         4.9         3.0          1.4         0.2  setosa
   3 │         4.7         3.2          1.3         0.2  setosa
   4 │         4.6         3.1          1.5         0.2  setosa
   5 │         5.0         3.6          1.4         0.2  setosa
   6 │         5.4         3.9          1.7         0.4  setosa
   7 │         4.6         3.4          1.4         0.3  setosa
   8 │         5.0         3.4          1.5         0.2  setosa
  ⋮  │      ⋮           ⋮            ⋮           ⋮           ⋮
 144 │         6.8         3.2          5.9         2.3  virginica
 145 │         6.7         3.3          5.7         2.5  virginica
 146 │         6.7         3.0          5.2         2.3  virginica
 147 │         6.3         2.5          5.0         1.9  virginica
 148 │         6.5         3.0          5.2         2.0  virginica
 149 │         6.2         3.4          5.4         2.3  virginica
 150 │         5.9         3.0          5.1         1.8  virginica
                                                   135 rows omitted
@df iris marginalhist(:PetalLength, :PetalWidth)

2.3 Lorenz

using DifferentialEquations, Plots
function lorenz(du,u,p,t)
 du[1] = p[1]*(u[2]-u[1])
 du[2] = u[1]*(p[2]-u[3]) - u[2]
 du[3] = u[1]*u[2] - p[3]*u[3]
end
lorenz (generic function with 1 method)

u0 = [1., 5., 10.]
3-element Vector{Float64}:
  1.0
  5.0
 10.0
tspan = (0., 100.)
(0.0, 100.0)
p = (10.0,28.0,8/3)
(10.0, 28.0, 2.6666666666666665)

prob = ODEProblem(lorenz, u0, tspan,p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: 3-element Vector{Float64}:
  1.0
  5.0
 10.0

sol = solve(prob)
retcode: Success
Interpolation: automatic order switching interpolation
t: 1345-element Vector{Float64}:
   0.0
   0.0354861341350177
   0.060663987304508726
   0.10188870843348657
   0.1448498679508511
   0.19835698662820245
   0.2504996228815297
   0.3056774805827286
   0.35452829390108126
   0.40770996764610945
   ⋮
  99.45173791784457
  99.518021374309
  99.59097978439264
  99.67312000769799
  99.76097062755619
  99.83799900237096
  99.91002674077228
  99.97332805621275
 100.0
u: 1345-element Vector{Vector{Float64}}:
 [1.0, 5.0, 10.0]
 [2.315652235826171, 5.897559436807754, 9.406792559102836]
 [3.2377953504336423, 7.041031570475975, 9.233678311348147]
 [4.993868184294182, 9.832941470623643, 9.626109614367385]
 [7.421184550586369, 13.9492707343288, 11.582332131961147]
 [11.45976330054411, 19.753113924253217, 18.104295519177246]
 [15.476108075765072, 21.510870676754823, 29.887267106518163]
 [16.44746410789155, 13.124038426772492, 40.97121918674213]
 [12.877766779212378, 2.618867895780566, 41.25247282418837]
 [7.13696497785409, -3.0934174660708695, 35.50512864260428]
 ⋮
 [-0.5025151204275446, -0.9233106458810417, 12.73214191014436]
 [-0.8622297993958143, -1.5678646795275535, 10.720273913657204]
 [-1.596977328543456, -2.9915709013611846, 9.007157621516578]
 [-3.3449989980630654, -6.426103674780766, 8.073170709902875]
 [-7.491146016463139, -14.158521227220465, 10.705626018644528]
 [-13.727954455531295, -22.207184489881456, 22.811039405817027]
 [-17.102797974189485, -15.688127772886453, 40.35905720289506]
 [-12.545627247564829, -1.2532166305907384, 41.513591364025444]
 [-9.430181940358784, 2.288447905782186, 38.488008826257456]

2.3.1 Static Plot

plot(sol, plotdensity=10000,lw=1.5)

plot(sol, plotdensity=10000, vars=(1,2))

plot(sol, plotdensity=10000, vars=(1,3))

plot(sol, plotdensity=10000, vars=(2,3))

plot(sol, plotdensity=10000, vars=(1,2,3))

2.3.2 Dynamic Plot

2.4 Wave Equation

using DifferentialEquations
struct Wave
    f
    g
    c
end

2.4.1 function plot_wave

using Plots

function plot_wave(y, file_name; fps = 60, kwargs...)
    anim = @animate for (i, y_row) in enumerate(eachrow(y))
        plot(
            y_row;
            title = "t = $(i-1)Δt",
            xlabel = "x",
            ylabel = "y(t, x)",
            legend = false,
            linewidth = 2,
            linecolor = "red",
            kwargs...
        )
    end
    gif(anim, file_name; fps, show_msg = false)

    return nothing
end
plot_wave (generic function with 1 method)

2.4.2 function solve_wave

function solve_wave(T, L, wave::Wave; n_t=100, n_x=100)
    ts = range(0, T; length=n_t)
    xs = range(0, L; length=n_x)
    Δt = ts[2] - ts[1]
    Δx = xs[2] - xs[1]
    y = zeros(n_t, n_x)

    # boundary conditions
    y[:,1] .= wave.f(0)
    y[:,end] .= wave.f(L)

    # initial conditions
    y[1,2:end-1] = wave.f.(xs[2:end-1])
    y[2,2:end-1] = y[1,2:end-1] + Δt*wave.g.(xs[2:end-1])

    # solution for t = 2Δt, 3Δt, ..., T
    for t in 2:n_t-1, x in 2:n_x-1
        ∂y_xx = (y[t, x+1] - 2*y[t, x] + y[t, x-1])/Δx^2
        y[t+1, x] = c^2 * Δt^2 * ∂y_xx  + 2*y[t, x] - y[t-1, x]
    end

    return y
end
solve_wave (generic function with 1 method)

2.4.3 Generate Wave \(f(x,L)\)

f(x,L) = 3*cos(x)*exp(-(x-L/2)^2) + x/L - sin(x+ L/T)
f (generic function with 1 method)
g(x) = 0
g (generic function with 1 method)

L = 2.0*pi
6.283185307179586
T = 500
500
c = 0.02
0.02

2.4.4 Display wave1.gif

wave = Wave(x -> f(x,L), g, c)
Wave(var"#61#62"(), g, 0.02)

y1 = solve_wave(T, L, wave; n_t=501, n_x=101)
501×101 Matrix{Float64}:
 -0.0124109  -0.0650979  -0.117457  …  1.09319   1.04047   0.987589
 -0.0124109  -0.0650979  -0.117457     1.09319   1.04047   0.987589
 -0.0124109  -0.0650646  -0.117397     1.09315   1.04046   0.987589
 -0.0124109  -0.0649988  -0.117278     1.09307   1.04042   0.987589
 -0.0124109  -0.0649017  -0.1171       1.09295   1.04037   0.987589
 -0.0124109  -0.0647751  -0.116863  …  1.09278   1.04029   0.987589
 -0.0124109  -0.0646204  -0.116567     1.09258   1.0402    0.987589
 -0.0124109  -0.0644392  -0.116214     1.09233   1.04008   0.987589
 -0.0124109  -0.0642322  -0.115804     1.09204   1.03994   0.987589
 -0.0124109  -0.064      -0.115339     1.09171   1.03977   0.987589
  ⋮                                 ⋱                      ⋮
 -0.0124109  -0.124486   -0.231656     0.497841  0.740035  0.987589
 -0.0124109  -0.129982   -0.242579     0.488933  0.735548  0.987589
 -0.0124109  -0.134974   -0.252518     0.481101  0.731611  0.987589
 -0.0124109  -0.139458   -0.261462  …  0.474354  0.728228  0.987589
 -0.0124109  -0.143431   -0.269406     0.468697  0.725401  0.987589
 -0.0124109  -0.146893   -0.276349     0.464129  0.723129  0.987589
 -0.0124109  -0.149845   -0.282295     0.460647  0.721411  0.987589
 -0.0124109  -0.152293   -0.287251     0.458241  0.720241  0.987589
 -0.0124109  -0.154241   -0.291229  …  0.456899  0.719613  0.987589
plot_wave(y1, "wave1.gif"; ylims=(-3,4))

2.4.5 Display wave2.gif

y2 = solve_wave(T, L, wave; n_t=501, n_x=200)
501×200 Matrix{Float64}:
 -0.0124109  -0.0389118  -0.065362   …  1.04074   1.01417   0.987589
 -0.0124109  -0.0389118  -0.065362      1.04074   1.01417   0.987589
 -0.0124109  -0.0388915  -0.0653286     1.04072   1.01417   0.987589
 -0.0124109  -0.0388537  -0.0652618     1.04069   1.01415   0.987589
 -0.0124109  -0.0388021  -0.0651628     1.04064   1.01413   0.987589
 -0.0124109  -0.0387382  -0.0650342  …  1.04056   1.01409   0.987589
 -0.0124109  -0.0386618  -0.0648792     1.04047   1.01404   0.987589
 -0.0124109  -0.0385719  -0.0646992     1.04034   1.01398   0.987589
 -0.0124109  -0.0384686  -0.0644935     1.0402    1.01391   0.987589
 -0.0124109  -0.038352   -0.0642611     1.04003   1.01382   0.987589
  ⋮                                  ⋱                      
 -0.0124109  -0.0696018  -0.126165      0.737736  0.86232   0.987589
 -0.0124109  -0.0721871  -0.13133       0.733589  0.860243  0.987589
 -0.0124109  -0.0745181  -0.135989      0.729993  0.858444  0.987589
 -0.0124109  -0.076594   -0.140141   …  0.72695   0.856923  0.987589
 -0.0124109  -0.0784145  -0.143784      0.724462  0.85568   0.987589
 -0.0124109  -0.0799805  -0.14692       0.722524  0.854714  0.987589
 -0.0124109  -0.0812936  -0.149554      0.721134  0.854023  0.987589
 -0.0124109  -0.082357   -0.151692      0.720283  0.853603  0.987589
 -0.0124109  -0.083175   -0.153341   …  0.719963  0.853451  0.987589
plot_wave(y2, "wave2.gif"; ylims=(-3,4) )

2.4.6 Display wave3.gif

y3 = solve_wave(T, L, wave; n_t=501, n_x=300)
501×300 Matrix{Float64}:
 -0.0124109  -0.0300526  -0.0476768  …  1.02297   1.00528   0.987589
 -0.0124109  -0.0300526  -0.0476768     1.02297   1.00528   0.987589
 -0.0124109  -0.0300367  -0.0476521     1.02296   1.00528   0.987589
 -0.0124109  -0.0300113  -0.0476028     1.02294   1.00527   0.987589
 -0.0124109  -0.0299779  -0.0475348     1.02291   1.00525   0.987589
 -0.0124109  -0.0299353  -0.0474503  …  1.02286   1.00523   0.987589
 -0.0124109  -0.0298842  -0.0473478     1.0228    1.0052    0.987589
 -0.0124109  -0.0298243  -0.047228      1.02271   1.00516   0.987589
 -0.0124109  -0.0297556  -0.0470907     1.02262   1.00511   0.987589
 -0.0124109  -0.0296782  -0.0469358     1.0225    1.00505   0.987589
  ⋮                                  ⋱                      
 -0.0124109  -0.0505706  -0.0885455     0.820852  0.90412   0.987589
 -0.0124109  -0.0522692  -0.091941      0.818135  0.90276   0.987589
 -0.0124109  -0.0537988  -0.0949991     0.815785  0.901585  0.987589
 -0.0124109  -0.0551584  -0.0977185  …  0.813805  0.900595  0.987589
 -0.0124109  -0.0563482  -0.100099      0.812193  0.89979   0.987589
 -0.0124109  -0.0573689  -0.102142      0.810949  0.899168  0.987589
 -0.0124109  -0.058222   -0.103851      0.810068  0.898729  0.987589
 -0.0124109  -0.0589098  -0.105229      0.809546  0.89847   0.987589
 -0.0124109  -0.0594349  -0.106284   …  0.809378  0.898388  0.987589
plot_wave(y3, "wave3.gif"; ylims=(-3,4) )