|
Moderator
Join Date: Nov 2001
Location: Netherlands
Posts: 6,364
|
@madshi,
Start with the polynomials
Code:
Y_I(x) = g3*x^3 + g2*x^2 + g1*x + g0, for x in [-3,-2]
Y_I(x) = f3*x^3 + f2*x^2 + f1*x + f0, for x in [-2,-1]
Y_I(x) = e3*x^3 + e2*x^2 + e1*x + e0, for x in [-1,0]
Y_I(x) = a3*x^3 + a2*x^2 + a1*x + a0, for x in [0,1]
Y_I(x) = b3*x^3 + b2*x^2 + b1*x + b0, for x in [1,2]
Y_I(x) = c3*x^3 + c2*x^2 + c1*x + c0, for x in [2,3]
Y_I(x) = d3*x^3 + d2*x^2 + d1*x + d0, for x in [3,4]
The boundary conditions should result in:
Code:
> eqn:={
> y0 = -27*g3 + 9*g2 - 3*g1 + g0,
> y1 = -8*g3 + 4*g2 - 2*g1 + g0,
> y1 = -8*f3 + 4*f2 - 2*f1 + f0,
> y2 = -f3 + f2 - f1 + f0,
> y2 = -e3 + e2 - e1 + e0,
> y3 = e0,
> y3 = a0,
> y4 = a3 + a2 + a1 + a0,
> y4 = b3 + b2 + b1 + b0,
> y5 = 8*b3 + 4*b2 + 2*b1 + b0,
> y5 = 8*c3 + 4*c2 + 2*c1 + c0,
> y6 = 27*c3 + 9*c2 + 3*c1 + c0,
> y6 = 27*d3 + 9*d2 + 3*d1 + d0,
> y7 = 64*d3 + 16*d2 + 4*d1 + d0,
> 12*g3 - 4*g2 + g1 = 12*f3 - 4*f2 + f1,
> 3*f3 - 2*f2 + f1 = 3*e3 - 2*e2 + e1,
> e1 = a1,
> 3*a3 + 2*a2 + a1 = 3*b3 + 2*b2 + b1,
> 12*b3 + 4*b2 + b1 = 12*c3 + 4*c2 + c1,
> 27*c3 + 6*c2 + c1 = 27*d3 + 6*d2 + d1,
> -12*g3 + 2*g2 = -12*f3 + 2*f2,
> -6*f3 + 2*f2 = -6*e3 + 2*e2,
> 2*e2 = 2*a2,
> 6*a3 + 2*a2 = 6*b3 + 2*b2,
> 12*b3 + 2*b2 = 12*c3 + 2*c2,
> 18*c3 + 2*c2 = 18*d3 + 2*d2,
> -18*g3 + 2*g2 = 0,
> 24*d3 + 2*d2 = 0};
The remaining Maple calculations for Spline64:
Code:
> solution:=solve(eqn,{a0,a1,a2,a3,b0,b1,b2,b3,c0,c1,c2,c3,d0,d1,d2,d3,e0,e1,e2,e3,f0,f1,f2,f3,g0,g1,g2,g3});
> a0s:=subs(solution, a0); a1s:=subs(solution, a1); a2s:=subs(solution, a2); a3s:=subs(solution, a3);
a0s := y3
2340 156 624 26 97 582
a1s := ---- y4 + ---- y6 - ---- y5 - ---- y7 - ---- y0 + ---- y1
2911 2911 2911 2911 2911 2911
2328
- ---- y2 - 3/2911 y3
2911
4050 270 1080 45 168 1008
a2s := ---- y4 + ---- y6 - ---- y5 - ---- y7 + ---- y0 - ---- y1
2911 2911 2911 2911 2911 2911
4032 6387
+ ---- y2 - ---- y3
2911 2911
49 24
a3s := - -- y4 - 6/41 y6 + -- y5 + 1/41 y7 - 1/41 y0 + 6/41 y1
41 41
24 49
- -- y2 + -- y3
41 41
> w:= a3s*x^3 + a2s*x^2 + a1s*x + a0s;
/ 49 24
w := |- -- y4 - 6/41 y6 + -- y5 + 1/41 y7 - 1/41 y0 + 6/41 y1
\ 41 41
24 49 \ 3 /4050 270 1080 45
- -- y2 + -- y3| x + |---- y4 + ---- y6 - ---- y5 - ---- y7
41 41 / \2911 2911 2911 2911
168 1008 4032 6387 \ 2 /2340
+ ---- y0 - ---- y1 + ---- y2 - ---- y3| x + |---- y4
2911 2911 2911 2911 / \2911
156 624 26 97 582 2328
+ ---- y6 - ---- y5 - ---- y7 - ---- y0 + ---- y1 - ---- y2
2911 2911 2911 2911 2911 2911
\
- 3/2911 y3| x + y3
/
> w0:=coeff(w,y0); w1:=coeff(w,y1); w2:=coeff(w,y2); w3:=coeff(w,y3);
3 168 2 97
w0 := - 1/41 x + ---- x - ---- x
2911 2911
3 1008 2 582
w1 := 6/41 x - ---- x + ---- x
2911 2911
24 3 4032 2 2328
w2 := - -- x + ---- x - ---- x
41 2911 2911
49 3 6387 2
w3 := -- x - ---- x - 3/2911 x + 1
41 2911
> w4:=coeff(w,y4); w5:=coeff(w,y5); w6:=coeff(w,y6); w7:=coeff(w,y7);
49 3 4050 2 2340
w4 := - -- x + ---- x + ---- x
41 2911 2911
24 3 1080 2 624
w5 := -- x - ---- x - ---- x
41 2911 2911
3 270 2 156
w6 := - 6/41 x + ---- x + ---- x
2911 2911
3 45 2 26
w7 := 1/41 x - ---- x - ---- x
2911 2911
> w0+w1+w2+w3+w4+w5+w6+w7;
1
Note that (w4,w5,w6,w7) are mirrored versions of (w3,w2,w1,w0). You can see that by replacing (y0,y1,y2,y3) by (y7,y6,y5,y4) and x by (1-x).
Last edited by Wilbert; 16th May 2009 at 23:02.
|