381 - Clamped Spline using equations

379 days ago by Professor381

This is an example of a clamped spline where all of the equations are listed out in the resulting system of equations.  Uses exactly four points plus prescribed slopes at the two ends.

This data is from "TeaTime" question #9 in section 5.2.
xs = [0.1 0.2 0.3 0.4] ys = [-0.62 -0.28 0.0066 0.24] mleft = 0.5 mright = 0.1 plot(xs,ys,'+') global a1 b1 c1 d1 a2 b2 c2 d2 a3 b3 c3 d3 
        
xs =

 0.1 0.2 0.3 0.4

ys =

 -0.62 -0.28 0.0066 0.24

mleft = 0.5
mright = 0.1
                                                                        

  0.4
|--------------------------------------------------------------------|  

      |          +           +          +           +          +        
|   
      |                                                                 
|   
  0.2 |-+                                                               
+-|   
      |                                                                 
|   
      |                                                                 
|   
      |                                                                 
|   
    0 |-+                                           A                   
+-|   
      |                                                                 
|   
      |                                                                 
|   
 -0.2 |-+                                                               
+-|   
      |                      A                                          
|   
      |                                                                 
|   
 -0.4 |-+                                                               
+-|   
      |                                                                 
|   
      |                                                                 
|   
      |                                                                 
|   
 -0.6 |-+                                                               
+-|   
      |                                                                 
|   
      |          +           +          +           +          +        
|   
 -0.8
|--------------------------------------------------------------------|  

     0.1       0.15         0.2       0.25         0.3       0.35       
0.4  
                                                                        

xs =

 0.1 0.2 0.3 0.4

ys =

 -0.62 -0.28 0.0066 0.24

mleft = 0.5
mright = 0.1
                                                                               
  0.4 |--------------------------------------------------------------------|   
      |          +           +          +           +          +           |   
      |                                                                    |   
  0.2 |-+                                                                +-|   
      |                                                                    |   
      |                                                                    |   
      |                                                                    |   
    0 |-+                                           A                    +-|   
      |                                                                    |   
      |                                                                    |   
 -0.2 |-+                                                                +-|   
      |                      A                                             |   
      |                                                                    |   
 -0.4 |-+                                                                +-|   
      |                                                                    |   
      |                                                                    |   
      |                                                                    |   
 -0.6 |-+                                                                +-|   
      |                                                                    |   
      |          +           +          +           +          +           |   
 -0.8 |--------------------------------------------------------------------|   
     0.1       0.15         0.2       0.25         0.3       0.35         0.4
A = [1 xs(1) xs(1)^2 xs(1)^3 0 0 0 0 0 0 0 0 ys(1); 1 xs(2) xs(2)^2 xs(2)^3 0 0 0 0 0 0 0 0 ys(2); 0 0 0 0 1 xs(2) xs(2)^2 xs(2)^3 0 0 0 0 ys(2); 0 0 0 0 1 xs(3) xs(3)^2 xs(3)^3 0 0 0 0 ys(3); 0 0 0 0 0 0 0 0 1 xs(3) xs(3)^2 xs(3)^3 ys(3); 0 0 0 0 0 0 0 0 1 xs(4) xs(4)^2 xs(4)^3 ys(4); 0 1 2*xs(2) 3*xs(2)^2 0 -1 -2*xs(2) -3*xs(2)^2 0 0 0 0 0; 0 0 0 0 0 1 2*xs(3) 3*xs(3)^2 0 -1 -2*xs(3) -3*xs(3)^2 0; 0 0 2 6*xs(2) 0 0 -2 -6*xs(2) 0 0 0 0 0; 0 0 0 0 0 0 2 6*xs(3) 0 0 -2 -6*xs(3) 0; 0 1 2*xs(1) 3*xs(1)^2 0 0 0 0 0 0 0 0 mleft; 0 0 0 0 0 0 0 0 0 1 2*xs(4) 3*xs(4)^2 mright ] c = rref(A) a1 = c(1,13) b1 = c(2,13) c1 = c(3,13) d1 = c(4,13) a2 = c(5,13) b2 = c(6,13) c2 = c(7,13) d2 = c(8,13) a3 = c(9,13) b3 = c(10,13) c3 = c(11,13) d3 = c(12,13) 
        
A =

 1 0.1 0.01 0.001 0 0 0 0 0 0 0 0 -0.62
 1 0.2 0.04 0.008 0 0 0 0 0 0 0 0 -0.28
 0 0 0 0 1 0.2 0.04 0.008 0 0 0 0 -0.28
 0 0 0 0 1 0.3 0.09 0.027 0 0 0 0 0.0066
 0 0 0 0 0 0 0 0 1 0.3 0.09 0.027 0.0066
 0 0 0 0 0 0 0 0 1 0.4 0.16 0.064 0.24
 0 1 0.4 0.12 0 -1 -0.4 -0.12 0 0 0 0 0
 0 0 0 0 0 1 0.6 0.27 0 -1 -0.6 -0.27 0
 0 0 2 1.2 0 0 -2 -1.2 0 0 0 0 0
 0 0 0 0 0 0 2 1.8 0 0 -2 -1.8 0
 0 1 0.2 0.03 0 0 0 0 0 0 0 0 0.5
 0 0 0 0 0 0 0 0 0 1 0.8 0.48 0.1

c =

 1 0 0 0 0 0 0 0 0 0 0 0 0.110773
 0 1 0 0 0 0 0 0 0 0 0 0 -17.5693
 0 0 1 0 0 0 0 0 0 0 0 0 127.155
 0 0 0 1 0 0 0 0 0 0 0 0 -245.387
 0 0 0 0 1 0 0 0 0 0 0 0 -2.6744
 0 0 0 0 0 1 0 0 0 0 0 0 24.2083
 0 0 0 0 0 0 1 0 0 0 0 0 -81.7333
 0 0 0 0 0 0 0 1 0 0 0 0 102.76
 0 0 0 0 0 0 0 0 1 0 0 0 4.56736
 0 0 0 0 0 0 0 0 0 1 0 0 -48.2093
 0 0 0 0 0 0 0 0 0 0 1 0 159.659
 0 0 0 0 0 0 0 0 0 0 0 1 -165.453

a1 = 0.110773
b1 = -17.5693
c1 = 127.155
d1 = -245.387
a2 = -2.6744
b2 = 24.2083
c2 = -81.7333
d2 = 102.76
a3 = 4.56736
b3 = -48.2093
c3 = 159.659
d3 = -165.453
A =

 1 0.1 0.01 0.001 0 0 0 0 0 0 0 0 -0.62
 1 0.2 0.04 0.008 0 0 0 0 0 0 0 0 -0.28
 0 0 0 0 1 0.2 0.04 0.008 0 0 0 0 -0.28
 0 0 0 0 1 0.3 0.09 0.027 0 0 0 0 0.0066
 0 0 0 0 0 0 0 0 1 0.3 0.09 0.027 0.0066
 0 0 0 0 0 0 0 0 1 0.4 0.16 0.064 0.24
 0 1 0.4 0.12 0 -1 -0.4 -0.12 0 0 0 0 0
 0 0 0 0 0 1 0.6 0.27 0 -1 -0.6 -0.27 0
 0 0 2 1.2 0 0 -2 -1.2 0 0 0 0 0
 0 0 0 0 0 0 2 1.8 0 0 -2 -1.8 0
 0 1 0.2 0.03 0 0 0 0 0 0 0 0 0.5
 0 0 0 0 0 0 0 0 0 1 0.8 0.48 0.1

c =

 1 0 0 0 0 0 0 0 0 0 0 0 0.110773
 0 1 0 0 0 0 0 0 0 0 0 0 -17.5693
 0 0 1 0 0 0 0 0 0 0 0 0 127.155
 0 0 0 1 0 0 0 0 0 0 0 0 -245.387
 0 0 0 0 1 0 0 0 0 0 0 0 -2.6744
 0 0 0 0 0 1 0 0 0 0 0 0 24.2083
 0 0 0 0 0 0 1 0 0 0 0 0 -81.7333
 0 0 0 0 0 0 0 1 0 0 0 0 102.76
 0 0 0 0 0 0 0 0 1 0 0 0 4.56736
 0 0 0 0 0 0 0 0 0 1 0 0 -48.2093
 0 0 0 0 0 0 0 0 0 0 1 0 159.659
 0 0 0 0 0 0 0 0 0 0 0 1 -165.453

a1 = 0.110773
b1 = -17.5693
c1 = 127.155
d1 = -245.387
a2 = -2.6744
b2 = 24.2083
c2 = -81.7333
d2 = 102.76
a3 = 4.56736
b3 = -48.2093
c3 = 159.659
d3 = -165.453
function y = S1(x,a1,b1,c1,d1) y = a1 + b1 .*x + c1*x .^2 + d1*x .^3 end function y = S2(x,a2,b2,c2,d2) y = a2 + b2 .*x + c2*x .^2 + d2*x .^3 end function y = S3(x,a3,b3,c3,d3) y = a3 + b3 .*x + c3*x .^2 + d3*x .^3 end 
       
S1(0.1,a1,b1,c1,d1) 
        
y = -0.62
ans = -0.62
y = -0.62
ans = -0.62
S2(0.3,a2,b2,c2,d2) 
        
y = 0.0066
ans = 0.0066
y = 0.0066
ans = 0.0066
S3(0.4,a3,b3,c3,d3) 
        
y = 0.24
ans = 0.24
y = 0.24
ans = 0.24

Now, let's put this all together into one cell that we can paste into sagecell.sagemath.org......
xs = [0.1 0.2 0.3 0.4] ys = [-0.62 -0.28 0.0066 0.24] plot(xs,ys,'+') hold on global a1 b1 c1 d1 a2 b2 c2 d2 a3 b3 c3 d3 A = [1 xs(1) xs(1)^2 xs(1)^3 0 0 0 0 0 0 0 0 ys(1); 1 xs(2) xs(2)^2 xs(2)^3 0 0 0 0 0 0 0 0 ys(2); 0 0 0 0 1 xs(2) xs(2)^2 xs(2)^3 0 0 0 0 ys(2); 0 0 0 0 1 xs(3) xs(3)^2 xs(3)^3 0 0 0 0 ys(3); 0 0 0 0 0 0 0 0 1 xs(3) xs(3)^2 xs(3)^3 ys(3); 0 0 0 0 0 0 0 0 1 xs(4) xs(4)^2 xs(4)^3 ys(4); 0 1 2*xs(2) 3*xs(2)^2 0 -1 -2*xs(2) -3*xs(2)^2 0 0 0 0 0; 0 0 0 0 0 1 2*xs(3) 3*xs(3)^2 0 -1 -2*xs(3) -3*xs(3)^2 0; 0 0 2 6*xs(2) 0 0 -2 -6*xs(2) 0 0 0 0 0; 0 0 0 0 0 0 2 6*xs(3) 0 0 -2 -6*xs(3) 0; 0 1 2*xs(1) 3*xs(1)^2 0 0 0 0 0 0 0 0 0.5; 0 0 0 0 0 0 0 0 0 1 2*xs(4) 3*xs(4)^2 0.1 ] c = rref(A) a1 = c(1,13) b1 = c(2,13) c1 = c(3,13) d1 = c(4,13) a2 = c(5,13) b2 = c(6,13) c2 = c(7,13) d2 = c(8,13) a3 = c(9,13) b3 = c(10,13) c3 = c(11,13) d3 = c(12,13) function y = S1(x,a1,b1,c1,d1) y = a1 + b1 .*x + c1*x .^2 + d1*x .^3 end function y = S2(x,a2,b2,c2,d2) y = a2 + b2 .*x + c2*x .^2 + d2*x .^3 end function y = S3(x,a3,b3,c3,d3) y = a3 + b3 .*x + c3*x .^2 + d3*x .^3 end u1 = 0.1:0.005:0.2; v1 = S1(u1,a1,b1,c1,d1) plot(u1,v1) u2 = 0.2:0.005:0.3; v2 = S2(u2,a2,b2,c2,d2) plot(u2,v2) u3 = 0.3:0.005:0.4; v3 = S3(u3,a3,b3,c3,d3) plot(u3,v3)