In [1]:

    

using Symata, NMfE, Plots


    

In [2]:

    

isymata()


    

    
Out[2]:





$$ isymata \!  \left(  \right)  $$

In [3]:

    

ClearAll(xi, yi, Y, R, sol, ytilde)
xi = [0, 1//2, 1]
yi = [0, a, 1]
Y(x_) := LagrangePolynomial(xi, yi)
#
# Can be formulated as ytile(x) = F(x) + C1(a) * Ψ(x)
#
# F(x_) := 2*x^2 - x
# C(a_) := -4a
# Ψ(x_) := x^2 - x
#
#Y(x_) := F(x) + C(a)*Ψ(x)
#
R(x_) = Simplify(D(Y(x), x, 2) - 3*x - 4*Y(x))
R(x_) = Simplify(R(x) ./ (a => -C1/4))
SetJ(r, ToString(Simplify(R(x))))


    

    
Out[3]:





$$ \text{"4.0 + 2.0C1 - 3x - 4x*(-1.0 + 2.0x + C1*(-1 + x))"} $$

In [4]:

    

Julia();


    

In [5]:

    

isymata()


    

In [6]:

    

Println("\nY(x) = ", Y(x), "\n")
Println("R(x) = ", R(x), "\n")


    

    




Y(x) = x*(-1.0 + 2.0x - 4.0a*(-1 + x))

R(x) = 4.0 + 2.0C1 - 3x - 4x*(-1.0 + 2.0x + C1*(-1 + x))

In [7]:

    

Julia()


    

    
Out[7]:




true

In [8]:

    

@assert r == "4.0 + 2.0C1 - 3x - 4x*(-1.0 + 2.0x + C1*(-1 + x))"


    

In [9]:

    

y_exact(x) = 7/4*(exp(2x)-exp(-2x))/(exp(2)-exp(-2))-3/4*x


    

    
Out[9]:





y_exact (generic function with 1 method)

In [26]:

    

ytildes_1pt = [
    (x) -> ytilde_1pt_collocation(x) - y_exact(x),
    (x) -> ytilde_1pt_subdomain(x) - y_exact(x),
    (x) -> ytilde_1pt_leastsquares(x) - y_exact(x),
    (x) -> ytilde_1pt_galerkin(x) - y_exact(x)
];

ytildes_2pt = [
    (x) -> ytilde_2pt_collocation(x) - y_exact(x),
    (x) -> ytilde_2pt_subdomain(x) - y_exact(x),
    (x) -> ytilde_2pt_leastsquares(x) - y_exact(x),
    (x) -> ytilde_2pt_galerkin(x) - y_exact(x)
];


    

    
Out[26]:





4-element Array{Function,1}:
 #13
 #14
 #15
 #16

In [ ]:

    

labels_1pt = ["ytilde_1pt_collocation", "ytilde_1pt_subdomain",
"ytilde_1pt_leastsquares", "ytilde_1pt_galerkin"];


    

In [ ]:

    

p1 = plot()
for (i, f) in enumerate(ytildes_1pt)
    lab = labels_1pt[i]
    if i == 1
      p1 = plot(f, 0, 1, label=lab, ylabel="ÿ - y")
    else
      plot!(p1, f, 0, 1, label=lab)
    end
end


    

In [ ]:

    

labels_2pt = ["ytilde_2pt_collocation", "ytilde_2pt_subdomain",
"ytilde_2pt_leastsquares", "ytilde_2pt_galerkin"];


    

In [ ]:

    

p2 = plot()
for (i, f) in enumerate(ytildes_2pt)
    lab = labels_2pt[i]
    if i == 1
      p2 = plot(f, 0, 1, label=lab, xlabel="x", ylabel="ÿ - y")
    else
      plot!(p2, f, 0, 1, label=lab)
    end
end


    

In [ ]:

    

plot(p1, p2, layout=(2, 1), link=:both, leg=true)