Example
Example
In[78]:=
u={0,.1125,.2763,.3794}
h={.1,.15,.1}
Out[78]=
{0, 0.1125, 0.2763, 0.3794}
Out[79]=
{0.1, 0.15, 0.1}
From our previous efforts (make sure you see where these came from!) , our
equations are:
.2s0+.1s1=3(.1125)=.3375
10s0 + 100/3 s1 + 20/3 s2
= 3/(.15)^2 (.2763) + (3/(.15)^2 - 3/(.1)^2)(.1125) - 3/(.1)^2(0)
=18.09
1/(.15)^2 s1 + (2/(.15)^2 + 2/(.1)^2)s2 + 1/(.1)^2 s3
= 3/(.1)^2 (.3794) + (3/(.15)^2 - 3/(.1)^2)(.2763) - 3/(.15)^2 (.1125)
=52.77
.1s2+.2s3=3(.1031)=.3093
In[80]:=
matrix={{2/10,1/10,0,0},{10,100/3,20/3,0},
{0,(20/3)^2,2(20/3)^2-300,100},{0,0,1/10,2/10}};
MatrixForm[matrix]
Out[80]=
1 1
- --
5 10 0 0
100 20
--- --
10 3 3 0
400 1900
--- -(----)
0 9 9 100
1 1
-- -
0 0 10 5
In[81]:=
du={.3375,18.09,52.77,.3039}
Out[81]=
{0.3375, 18.09, 52.77, 0.3039}
In[82]:=
s=LinearSolve[matrix,du]
Out[82]=
{1.70985, -0.0447011, 0.37223, 1.33339}
This gives s0, s1, s2, s3, which Mathematica sees as s[[1]], s[[2]], s[[3]], and s[[4]]
To the left of (0,0) we have a line with slope s0=1.70985 through (0,0).
In[83]:=
Clear[x]
In[84]:=
newspline0[x_]:=s[[1]]x;
newspline0[x]
Out[84]=
1.70985 x
Between x=0 and x=0.1 we have
In[85]:=
newspline1[x_]:=genherm[x,0,s[[1]],.1125,s[[2]],h[[1]]];
newspline1[x]
Out[85]=
-18 2 3
1.70985 x - 2.71051 10 x - 58.4851 x
Between x=.1 and x=.25 we have
In[86]:=
newspline2[x_]:=gc2[x,.1125,s[[2]],.2763,s[[3]]]/.{h1->h[[1]],h2->h[[2]]};
newspline2[x]
Out[86]=
2 3
0.399025 - 6.51089 x + 44.7074 x - 82.5098 x
In our Theory section we did not go beyond gc2, so we must find gc3 (spline3) to be between x=.25 and x=.35.
In[87]:=
Clear[u]
DSolve[{u''''[x]==0,u[.25]==.2763,u'[.25]==s[[3]],
u[.35]==.3794,u'[.35]==s[[4]]},u[x],x]/.x->x
Out[87]=
2 3
{{u[x] -> 1.37457 - 11.3858 x + 36.8804 x - 35.6385 x }}
In[88]:=
newspline3[x_]:=
1.374566719314079424 - 11.38577084837545128*x +
36.88044584837545133*x^2 - 35.63851985559566794*x^3
To the right of x=.35 we have a line with slope 1.33339 through (.35, .3795)
In[89]:=
newspline4[x_]:=.3795+s[[4]](x-.35);
newspline4[x]
Out[89]=
0.3795 + 1.33339 (-0.35 + x)
Let's put it together!
In[90]:=
newspline[x_]:=If[x<=0,newspline0[x],
If[x<=.1,newspline1[x],
If[x<=.25,newspline2[x],
If[x<=.35,newspline3[x],newspline4[x]]
]
]
]
In[91]:=
xvals={0,.1,.25,.35}
u={0,.1125,.2763,.3794}
Out[91]=
{0, 0.1, 0.25, 0.35}
Out[92]=
{0, 0.1125, 0.2763, 0.3794}
In[93]:=
plot3=ListPlot[Table[{xvals[[i]],u[[i]]},{i,1,4}],PlotStyle->{RGBColor[1,0,0],PointSize[.02]}];
plot4=Plot[newspline[x],{x,-.1,.4}];
Show[plot3,plot4];
The points that we have been considering actually are approximations from the
error function Erf[x] which is defined by
Erf[x]=2/Sqrt[Pi] Integral from 0 to x of Exp[-t^2]dt.
Let's see how we did.
In[94]:=
plot5=Plot[Erf[x],{x,-.1,.4}];
Show[plot3,plot4,plot5];
