-
Notifications
You must be signed in to change notification settings - Fork 3
/
Ex.7.13.function.jl
78 lines (65 loc) · 2.26 KB
/
Ex.7.13.function.jl
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
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
using NumericalMethodsforEngineers
using Test
@sym begin
CollocationWeightedResidualMethod(xi_, yi_, doprint_) := Module([Y, R],
begin
Y(x_) := LagrangePolynomial(xi, yi)
If(doprint, Println("\nY(x) = ", Y(x), "\n"))
#
# 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))
If(doprint, Println("R(x) = ", R(x), "\n"))
Return(R(x))
end
)
end
function wrm(x1,y1, doprint=true)
setsymata(:xin, List(x1...))
setsymata(:yin, List(y1...))
setsymata(:doprint, doprint ? true : false)
symparseeval("Compile(Evaluate(CollocationWeightedResidualMethod(xin,yin,doprint)))")
end
println("\nExample 7.13: y'' = 3x + 4y, y(0)=0, y(1)=1")
println("Residual for Weighted Residual Method using 1 point Lagragian Polynomial\n")
#
# A consequence of the Julia function based approach is that Y(x) and R(x) are local
# to the Symata function and are not available for examination, e.g. @sym R(x)
#
lp_eq = @sym ToString(LagrangePolynomial([0.0,0.5,1.0], [0.0,a,1.0]))
println(lp_eq)
@test lp_eq == "x*(-1.0 + 2.0*x + (-4.0)*a*(-1.0 + x))"
wrm_eq = @sym ToString(CollocationWeightedResidualMethod([0.0,0.5,1.0], [0.0,a,1.0], false))
println(wrm_eq)
@test wrm_eq == "4.0 + 2.0*C1 - 3x - 4x*(-1.0 + 2.0*x + C1*(-1.0 + x))"
#
# Another consideration is that below yt is yt(C1::Any, x::Any)
#
yt = wrm([0.0,0.5,1.0], [0.0,:a,1.0], true)
println("( Example 7.13 gives: R = -4x^2*(2 + C1) + x*(1 + 4C1) + 2*(2 + C1) )", "\n")
"
```
The function ytilde(x, C1) contains the 1-point Weighted Residual Method
polynomial with parameter C1.
```
Arguments to ytilde(x, C1):
x::Float64 : x value\n
C1::Float64 : Weighted Residual Method constant.
C1 can be estimated by e.g. collocation, subdomain, least squares or Galerkin.
See examples 7.15 through 7.18 and Fig.7.14.jl
"
ytilde(x, C1) = yt(C1, x)
println("ytilde(1.0, 0.0) = $(ytilde(1.0, 0.0))")
@test ytilde(1.0, 0.0) == -3.0
println("ytilde(0.0, 1.0) = $(ytilde(0.0, 1.0))")
@test ytilde(0.0, 1.0) == 6.0
#
# For now I will continue to use the E.7.13.jl approach.
#