/*WORK:"202107-1MaximaTrial.html";
: 　　トピックス１
　　　. 関数列と極限関数
◆コマンド番号202107-11
●取得データの、差分商による補間（ニュートン）
●●●●●●●●●●●●●●●●*/
/*(kill(all), print("インストールMaximaでは、先頭の／*を、次行の先頭にカット＆ペーストすること"))$
*/

(MONDAI:"・Fuse202107-1　(1,1),(2,½),(4,¼),(7,1/7)をグラフが通る関数を、3次関数で近似し、y=1/x との誤差を検討せよ。"/*(◆修正部分１)*/,
KOUSATU:"ニュートンの補間法で考えると、・・・"/**/)$

(KATEI:"SetKATEI",
subl:[x1=1,x2=2,x3=4,x4=7,f1=1,f2=1/2,f3=1/4,f4=1/7],print(subl)/*(◆修正部分２)*/,
KATEI:"EndKATEI")$

(Kansuu0:"",
I(x,a):=if(x> a) then 1 else 0,
gr(T,c,l,pP,pK,lab,xm,xM,ym,yM):=block(
/**/wxdraw2d(
/*....*/title =T, yrange=[ym,yM] /* "タイトルと値域を確認し" */,
/*....*/color=blue,line_width=1,
/*....*/implicit(c,x,xm,xM,y,ym,yM),
/*....*/color=cyan,line_width=1,
/*....*/implicit(l,x,xm,xM,y,ym,yM),
/*....*/color=red,point_type =7,line_width=1,
/*....*/points_joined =true,point_size =1,points(pK[2]) /* "点pKを赤で打って、赤線で結んで" */,
/*....*/color=cyan,point_type =7,line_width=1,
/*....*/points_joined =true,point_size =1,points(pK[1]) /* "点pKをシアンで打って、シアン線で結んで" */,
/*....*/color=black,point_type =1,line_width=1,
/*....*/points_joined =false,point_size =2,points(pP) /* "点pPを青で打って、青線で結んで" */,
/*....*/color=red,lab[2][1],lab[2][2],lab[2][3],lab[2][4],lab[2][5],lab[2][6],lab[2][7],lab[2][8],lab[2][9],lab[2][10],lab[2][11],lab[2][12],lab[2][13],
/*....*/color=black,line_width=1,
/*....*/implicit(y=0,x,xm,xM,y,ym,yM),implicit(x=0,x,xm,xM,y,ym,yM) /* "x,y軸を引く" */,
/*....*/color=black)),
lab(nam,p):=label([nam,p[1],p[2]]) /* 図中の文字 */,
Kansuu0:"EndKansuu")$

(KANSUU:"Teigi",
X:[x1,x2,x3,x4],print("X =",subst(subl,X)),
f:[f1,f2,f3,f4],print("f =",subst(subl,f)),
d1:[(f[2]-f[1])/(X[2]-X[1]),(f[3]-f[2])/(X[3]-X[2]),(f[4]-f[3])/(X[4]-X[3])],print("d1 =",subst(subl,d1)),
d2:[(d1[2]-d1[1])/(X[3]-X[1]),(d1[3]-d1[2])/(X[4]-X[2])],print("d2 =",subst(subl,d2)),
d3:[(d2[2]-d2[1])/(X[4]-X[1])],print("d3 =",subst(subl,d3)),
f(x):=rat(subst(subl,(f[1] +d1[1]*(x-x1) +d2[1]*(x-x1)*(x-x2) +d3[1]*(x-x1)*(x-x2)*(x-x3)))),print("f(x) =",f(x)),
KANSUU:"EndTeigi")$

(KAKUNIN:"StartKAKUNIN　上記3次関数f(x)が、f(x1)=y1, f(x2)=y2,f(x3)=y3,f(x4)=y4 であることは、以下による。",
e10:["f(x1)=f[1]=f1","(定義・代入)",bfloat(subst(subl,[f(x1)=f[1],"=",f1]) )],print(e10),
e20:["f(x2)=f[1] +d1[1]*(x2-x1)","(定義・代入)",bfloat(subst(subl,[f(x2)=f[1] +d1[1]*(x2-x1)]) )],print(e20),
e21:[" =f1 +(f2-f1)/(x2-x1)*(x2-x1)","(代入)",bfloat(subst(subl,[" =",f1 +(f2-f1)/(x2-x1)*(x2-x1)]) )],print(e21),
e22:[" =f1 +(f2-f1) =f2","(簡約)",bfloat(subst(subl,[" =",f1 +(f2-f1)," =",f2]) )],print(e22),
e30:["f(x3)=f[1] +d1[1]*(x3-x1) +d2[1]*(x3-x1)*(x3-x2)","(定義・代入)",bfloat(subst(subl,[f(x3)=f[1] +d1[1]*(x3-x1) +d2[1]*(x3-x1)*(x3-x2)]) )],print(e30),
e31:[" =f1 +(f2-f1)/(x2-x1)*(x3-x1) +((f3-f2)/(x3-x2)-(f2-f1)/(x2-x1))/(x3-x1)*(x3-x1)*(x3-x2)","(代入)",
bfloat(subst(subl, [" =", f1 +(f2-f1)/(x2-x1)*(x3-x1) +((f3-f2)/(x3-x2)-(f2-f1)/(x2-x1))/(x3-x1)*(x3-x1)*(x3-x2)]) )],print(e31),
e32:[" =f1 +(f2-f1)/(x2-x1)*(x3-x1) +(f3-f2) -(f2-f1)/(x2-x1)*(x3-x2)","(簡約)",
bfloat(subst(subl, [" =",f1 +(f2-f1)/(x2-x1)*(x3-x1) +(f3-f2)-(f2-f1)/(x2-x1)*(x3-x2)]) )],print(e32),
e33:[" =f1 +(f2-f1)/(x2-x1)*((x3-x1) -(x3-x2)) +(f3-f2)","(簡約)",
bfloat(subst(subl, [" =",f1 +(f2-f1)/(x2-x1)*((x3-x1) -(x3-x2)) +(f3-f2)]) )],print(e33),
e34:[" =f1 +(f2-f1)/(x2-x1)*(x2-x1) +(f3-f2)","(簡約)", bfloat(subst(subl, [" =",f1 +(f2-f1)/(x2-x1)*(x2-x1) +(f3-f2)]) )],print(e34),
e35:[" =f1 +(f2-f1) +(f3-f2) =f3","(簡約)", bfloat(subst(subl, [" =",f1 +(f2-f1) +(f3-f2) =f3]) )],print(e35),
e40:["f(x4)=f[1] +d1[1]*(x4-x1) +d2[1]*(x4-x1)*(x4-x2) +d3[1]*(x4-x1)*(x4-x2)*(x4-x3)","(定義・代入)",
bfloat(subst(subl,[f(x4)=f[1] +d1[1]*(x4-x1) +d2[1]*(x4-x1)*(x4-x2) +d3[1]*(x4-x1)*(x4-x2)*(x4-x3)]) )],print(e40),
e41:[" =f1 +(f2-f1)/(x2-x1)*(x4-x1) +((f3-f2)/(x3-x2)-(f2-f1)/(x2-x1))/(x3-x1)*(x4-x1)*(x4-x2)
+(((f4-f3)/(x4-x3)-(f3-f2)/(x3-x2))/(x4-x2)-((f3-f2)/(x3-x2)-(f2-f1)/(x2-x1))/(x3-x1))/(x4-x1)*(x4-x1)*(x4-x2)*(x4-x3)","(代入)",
bfloat(subst(subl,[" =",f1 +(f2-f1)/(x2-x1)*(x4-x1) +((f3-f2)/(x3-x2)-(f2-f1)/(x2-x1))/(x3-x1)*(x4-x1)*(x4-x2)
+(((f4-f3)/(x4-x3)-(f3-f2)/(x3-x2))/(x4-x2)-((f3-f2)/(x3-x2)-(f2-f1)/(x2-x1))/(x3-x1))/(x4-x1)*(x4-x1)*(x4-x2)*(x4-x3)]) )],print(e41),
e42:[" =f1 +(f2-f1)/(x2-x1)*(x4-x1) +((f3-f2)/(x3-x2)-(f2-f1)/(x2-x1))/(x3-x1)*(x4-x1)*(x4-x2)
+(f4-f3)-(f3-f2)/(x3-x2)*(x4-x3)-((f3-f2)/(x3-x2)-(f2-f1)/(x2-x1))/(x3-x1)*(x4-x2)*(x4-x3)","(簡約)",
bfloat(subst(subl,[" =",f1 +(f2-f1)/(x2-x1)*(x4-x1) +((f3-f2)/(x3-x2)-(f2-f1)/(x2-x1))/(x3-x1)*(x4-x1)*(x4-x2)
+(f4-f3)-(f3-f2)/(x3-x2)*(x4-x3)-((f3-f2)/(x3-x2)-(f2-f1)/(x2-x1))/(x3-x1)*(x4-x2)*(x4-x3)]) )],print(e42),
e43:[" =f1 +(f2-f1)/(x2-x1)*((x4-x1)-(x4-x1)*(x4-x2)/(x3-x1)+(x4-x2)*(x4-x3)/(x3-x1))
+(f3-f2)/(x3-x2)*((x4-x1)*(x4-x2)/(x3-x1)-(x4-x3)-(x4-x2)*(x4-x3)/(x3-x1)) +(f4-f3)","(簡約)",
bfloat(subst(subl,[" =",f1 +(f2-f1)/(x2-x1)*((x4-x1)-(x4-x1)*(x4-x2)/(x3-x1)+(x4-x2)*(x4-x3)/(x3-x1))
+(f3-f2)/(x3-x2)*((x4-x1)*(x4-x2)/(x3-x1)-(x4-x3)-(x4-x2)*(x4-x3)/(x3-x1)) +(f4-f3)/**/]) )/**/],print(e43),
e44:[" =f1 +(f2-f1)/(x2-x1)*((x4-x1)-(x4-x2)/(x3-x1)*((x4-x1)-(x4-x3))
+(f3-f2)/(x3-x2)*((x4-x2)/(x3-x1)*((x4-x1)-(x4-x3)) -(x4-x3)) +(f4-f3)","(簡約)",
bfloat(subst(subl,[" =",f1 +(f2-f1)/(x2-x1)*((x4-x1)-(x4-x2)/(x3-x1)*((x4-x1)-(x4-x3)))
+(f3-f2)/(x3-x2)*((x4-x2)/(x3-x1)*((x4-x1)-(x4-x3)) -(x4-x3)) +(f4-f3)/**/]) )/**/],print(e44),
e45:[" =f1 +(f2-f1)/(x2-x1)*((x4-x1)-(x4-x2)) +(f3-f2)/(x3-x2)*((x4-x2)-(x4-x3)) +(f4-f3)","(簡約)",
bfloat(subst(subl,[" =",f1 +(f2-f1)/(x2-x1)*((x4-x1)-(x4-x2)) +(f3-f2)/(x3-x2)*((x4-x2)-(x4-x3)) +(f4-f3)/**/]) )/**/],print(e45),
e46:[" =f1 +(f2-f1) +(f3-f2) +(f4-f3) =f4","(簡約)",
bfloat(subst(subl,[" =",f1 +(f2-f1) +(f3-f2) +(f4-f3) =f4/**/]) )/**/],print(e46),
/**/END:"EndKAKUNIN")$
f(bfloat(300001));
sin((30+1/3600)/180*3.1415926535);

(GJUNBI:"グラフ準備",
xM:8, yM: 1.5, xm:-xM/8, ym:-yM/8, print([xm,xM,ym,yM]) /*出力範囲　(◆修正部分３)*/,

dA:[0,-0.1],dB:[0,-0.1],dC:[0,-0.1],dD:[0,-0.1] /**/,
dE:[-0.3,0.1],dF:[-0.3,0.1],dG:[-0.3,0.1],dH:[-0.3,0.1] /**/,
dO:[-0.3,-0.1],dP:[0.4,0],dQ:[0.4,0],dR:[0.4,0],dS:[0.4,0] /**/,

pA:[rhs(subl[1]),0],pB:[rhs(subl[2]),0],pC:[rhs(subl[3]),0],pD:[rhs(subl[4]),0] /**/,
pE:[0,rhs(subl[5])],pF:[0,rhs(subl[6])],pG:[0,rhs(subl[7])],pH:[0,rhs(subl[8])] /**/,
pP:[rhs(subl[1]),rhs(subl[5])],pQ:[rhs(subl[2]),rhs(subl[6])],pR:[rhs(subl[3]),rhs(subl[7])],pS:[rhs(subl[4]),rhs(subl[8])],pO:[0,0] /**/,

labA:lab("A",pA+dA),labB:lab("B",pB+dB),labC:lab("C",pC+dC),labD:lab("D",pD+dD) /**/,
labE:lab("E",pE+dE),labF:lab("F",pF+dF),labG:lab("G",pG+dG),labH:lab("H",pH+dH) /**/,
labP:lab("P",pP+dP),labQ:lab("Q",pQ+dQ),labR:lab("R",pR+dR),labS:lab("S",pS+dS),labO:lab("O",pO+dO) /**/,
pRB:[[pP,pA,pB,pQ,pB,pC,pR,pC,pD,pS],[pP,pE,pF,pQ,pF,pG,pR,pG,pH,pS]]/**/,print(float(pRB)) /*点出力準備*/,
lab:[13,[labA,labB,labC,labD,labE,labF,labG,labH,labP,labQ,labR,labS,labO]],print(float(lab)) /*ラベル準備*/,
pPG:[pO],print(float(pPG)) /**/,
/**/GJUNBI:"EndGJunbi")$

(GPrint:"グラフを書いて、視覚的に確認",
c:y =f(x) /**/,
T:"y =1/x (cyan)", c1: y=1/x /*(◆修正部分４)*/,
print([T,c,c1,pPG,pRB,lab,xm,xM,ym,yM]) /* 出力形式*/,
gr(T,c,c1,pPG,pRB,lab,xm,xM,ym,yM) /* グラフ出力すると*/,
/**/Gprint:"EndGPrint")$

/**//*◆◆◆◎上のMaximaコマンドを実行せよ。
◆◆◆◎上のコマンドを、下の4例（一番下は上と同じ）次により修正して実行し、結果を比較して考察せよ。
課題１　下の4つのうち、近似が最もうまくいっているのはどれか？
課題２　近似が、大きくずれる原因は何か？　
課題３　逆に、近似がうまくいく原因は何か？
課題４　通過する点を5つにして、4次関数で近似するにはどうするとよいか？
(KANSUU:"Teigi", ～　KANSUU:"EndTeigi")$　の部分を類推・試行錯誤してみよ。
例２
(MONDAI:"・202107-12　(1,1),(3,1/3),(5,1/5),(7,1/7)をグラフが通る関数を、3次関数で近似し、y=1/x との誤差を検討せよ。"/*(◆修正部分１)*/,
subl:[x1=1,x2=3,x3=5,x4=7,f1=1,f2=1/3,f3=1/5,f4=1/7],print(subl)/*(◆修正部分２)*/,
xM:8, yM: 1.5, xm:-xM/8, ym:-yM/8, print([xm,xM,ym,yM]) /*出力範囲　(◆修正部分３)*/,
T:"y =1/x (cyan)", c1:y =1/x /*(◆修正部分４)*/,

例３
(MONDAI:"・202107-13　(0,1),(1,2),(2,4),(3,8)をグラフが通る関数を、3次関数で近似し、y=2^x との誤差を検討せよ。"/*(◆修正部分１)*/,
subl:[x1=0,x2=1,x3=2,x4=3,f1=1,f2=2,f3=4,f4=8],print(subl)/*(◆修正部分２)*/,
xM:4, yM: 10, xm:-xM/8, ym:-yM/8, print([xm,xM,ym,yM]) /*出力範囲　(◆修正部分３)*/,
T:"y =2^x (cyan)", c1:y= 2^x /*(◆修正部分４)*/,

例４
(MONDAI:"・202107-14　(1,0),(2,1),(4,2),(8,3)をグラフが通る関数を、3次関数で近似し、y=log2(x) との誤差を検討せよ。"/*(◆修正部分１)*/,
subl:[x1=1,x2=2,x3=4,x4=8,f1=0,f2=1,f3=2,f4=3],print(subl)/*(◆修正部分２)*/,
xM:10, yM: 4, xm:-xM/8, ym:-yM/8, print([xm,xM,ym,yM]) /*出力範囲　(◆修正部分３)*/,
T:"y =log2(x) (cyan)", c1:y =log(x)/log(2)*I(x,0) /*(◆修正部分４)*/,

例５
(MONDAI:"・202107-15　(0,0),(10,sin(10/180*3.1415926535)),(20,sin(20/180*3.1415926535)),(30,0.5)をグラフが通る関数を、3次関数で近似し、y=sin(x) との誤差を検討せよ。"/*(◆修正部分１)*/,
subl:[x1=0,x2=10,x3=20,x4=30,f1=0,f2=sin(10/180*3.1415926535),f3=sin(20/180*3.1415926535),f4=0.5],print(subl)/*(◆修正部分２)*/,
xM:45, yM: 1, xm:-xM/8, ym:-yM/8, print([xm,xM,ym,yM]) /*出力範囲　(◆修正部分３)*/,
T:"y =sin(x) (cyan)", c1:y =sin(x/180*3.1415926535) /*(◆修正部分４)*/,

例１
(MONDAI:"・202107-11　(1,1),(2,½),(4,¼),(7,1/7)をグラフが通る関数を、3次関数で近似し、y=1/x との誤差を検討せよ。"/*(◆修正部分１)*/,
subl:[x1=1,x2=2,x3=4,x4=7,f1=1,f2=1/2,f3=1/4,f4=1/7],print(subl)/*(◆修正部分２)*/,
xM:8, yM: 1.5, xm:-xM/8, ym:-yM/8, print([xm,xM,ym,yM]) /*出力範囲　(◆修正部分３)*/,
T:"y =1/x (cyan)", c1: y=1/x /*(◆修正部分４)*/,

*/END:"Finish"$