Please I need Correction on this code particularly if I can make do without the declaration of vector in the third subroutine . The idea is to get maximum error. The code has 3 subroutine. The problem I think is in the third subroutine (Display of results).
Thank you in anticipation of positive response.
# First Declaration of the problem
restart:
Digits:=30:
interface(rtablesize=infinity):
f1:=proc(n)
y2[n]:
end proc:
f2:=proc(n)
-y1[n]+0.001*cos(t[n]):
end proc:
f3:=proc(n)
y4[n]:
end proc:
f4:=proc(n)
-y3[n]+0.001*sin(t[n]):
end proc:
F1:=proc(n)
f2(n):
end proc:
F2:=proc(n)
-(f1(n))-0.001*sin(t[n]):
end proc:
F3:=proc(n)
f4(n):
end proc:
F4:=proc(n)
-f3(n)+0.001*cos(t[n]):
end proc:
# Declaration of the Numerical methods
e1:=y1[n+2] = (7/23)*y1[n]+(16/23)*y1[n+1]+(12/23)*f1(n+2)*h+(16/23)*f1(n+1)*h-(2/23)*F1(n+2)*h^2+(2/23)*h*f1(n)+((24/3703)*y1[n]-(24/3703)*y1[n+1]+(48/18515)*f1(n+2)*h+(8/55545)*f1(n+1)*h-(116/55545)*F1(n+2)*h^2+(208/55545)*h*f1(n))*u^2+((901/2980915)*y1[n]-(901/2980915)*y1[n+1]+(7109/89427450)*f1(n+2)*h+(923/14904575)*f1(n+1)*h-(6241/89427450)*F1(n+2)*h^2+(14383/89427450)*h*f1(n))*u^4+((1979723/158376013950)*y1[n]-(1979723/158376013950)*y1[n+1]+(6364571/2375640209250)*f1(n+2)*h+(728327/215967291750)*f1(n+1)*h-(11785633/4751280418500)*F1(n+2)*h^2+(5106559/791880069750)*h*f1(n))*u^6+((6488435581/13259239887894000)*y1[n]-(6488435581/13259239887894000)*y1[n+1]+(8693517709/91794737685420000)*f1(n+2)*h+(260601208141/1789997384865690000)*f1(n+1)*h-(323357994149/3579994769731380000)*F1(n+2)*h^2+(891627999937/3579994769731380000)*h*f1(n))*u^8+((25090513463/1343541160668420000)*y1[n]-(25090513463/1343541160668420000)*y1[n+1]+(190450718149/55421072877572325000)*f1(n+2)*h+(47563947061/8210529315195900000)*f1(n+1)*h-(1475729910283/443368583020578600000)*F1(n+2)*h^2+(261738159769/27710536438786162500)*h*f1(n))*u^10+((244426606265778733/347060946154014557665200000)*y1[n]-(244426606265778733/347060946154014557665200000)*y1[n+1]+(1316372988977975777/10411828384620436729956000000)*f1(n+2)*h+(105391490263288387/473264926573656214998000000)*f1(n+1)*h-(1284959669761615073/10411828384620436729956000000)*F1(n+2)*h^2+(72506125749079249/204153497737655622156000000)*h*f1(n))*u^12:
e2:=h^2*F1(n+1) = (60/23)*y1[n]-(60/23)*y1[n+1]+(25/46)*f1(n+2)*h+(32/23)*f1(n+1)*h-(4/23)*F1(n+2)*h^2+(31/46)*h*f1(n)+((209/3703)*y1[n]-(209/3703)*y1[n+1]+(1313/222180)*f1(n+2)*h+(1304/55545)*f1(n+1)*h-(131/18515)*F1(n+2)*h^2+(6011/222180)*h*f1(n))*u^2+((77491/35770980)*y1[n]-(77491/35770980)*y1[n+1]+(574843/2146258800)*f1(n+2)*h+(113536/134141175)*f1(n+1)*h-(53461/178854900)*F1(n+2)*h^2+(2258041/2146258800)*h*f1(n))*u^4+((151508243/1900512167400)*y1[n]-(151508243/1900512167400)*y1[n+1]+(1290306599/114030730044000)*f1(n+2)*h+(18919693/647901875250)*f1(n+1)*h-(113769323/9502560837000)*F1(n+2)*h^2+(4470322013/114030730044000)*h*f1(n))*u^6+((42120775181/14464625332248000)*y1[n]-(42120775181/14464625332248000)*y1[n+1]+(332746636891/734357901483360000)*f1(n+2)*h+(302396120633/298332897477615000)*f1(n+1)*h-(369019384141/795554393273640000)*F1(n+2)*h^2+(13797329479621/9546652719283680000)*h*f1(n))*u^8+((18953368786273/177347433208231440000)*y1[n]-(18953368786273/177347433208231440000)*y1[n+1]+(2430202319484337/138330997902420523200000)*f1(n+2)*h+(310803544671199/8645687368901282700000)*f1(n+1)*h-(203453960588449/11527583158535043600000)*F1(n+2)*h^2+(7380568619069419/138330997902420523200000)*h*f1(n))*u^10+((16436168060905785763/4164731353848174691982400000)*y1[n]-(16436168060905785763/4164731353848174691982400000)*y1[n+1]+(167160345356705269819/249883881230890481518944000000)*f1(n+2)*h+(461636091223370027/354948694930242161248500000)*f1(n+1)*h-(13852288092290788813/20823656769240873459912000000)*F1(n+2)*h^2+(29059878239787610409/14699051837111204795232000000)*h*f1(n))*u^12:
e3:=y2[n+2] = (7/23)*y2[n]+(16/23)*y2[n+1]+(12/23)*f2(n+2)*h+(16/23)*f2(n+1)*h-(2/23)*F2(n+2)*h^2+(2/23)*h*f2(n)+((24/3703)*y2[n]-(24/3703)*y2[n+1]+(48/18515)*f2(n+2)*h+(8/55545)*f2(n+1)*h-(116/55545)*F2(n+2)*h^2+(208/55545)*h*f2(n))*u^2+((901/2980915)*y2[n]-(901/2980915)*y2[n+1]+(7109/89427450)*f2(n+2)*h+(923/14904575)*f2(n+1)*h-(6241/89427450)*F2(n+2)*h^2+(14383/89427450)*h*f2(n))*u^4+((1979723/158376013950)*y2[n]-(1979723/158376013950)*y2[n+1]+(6364571/2375640209250)*f2(n+2)*h+(728327/215967291750)*f2(n+1)*h-(11785633/4751280418500)*F2(n+2)*h^2+(5106559/791880069750)*h*f2(n))*u^6+((6488435581/13259239887894000)*y2[n]-(6488435581/13259239887894000)*y2[n+1]+(8693517709/91794737685420000)*f2(n+2)*h+(260601208141/1789997384865690000)*f2(n+1)*h-(323357994149/3579994769731380000)*F2(n+2)*h^2+(891627999937/3579994769731380000)*h*f2(n))*u^8+((25090513463/1343541160668420000)*y2[n]-(25090513463/1343541160668420000)*y2[n+1]+(190450718149/55421072877572325000)*f2(n+2)*h+(47563947061/8210529315195900000)*f2(n+1)*h-(1475729910283/443368583020578600000)*F2(n+2)*h^2+(261738159769/27710536438786162500)*h*f2(n))*u^10+((244426606265778733/347060946154014557665200000)*y2[n]-(244426606265778733/347060946154014557665200000)*y2[n+1]+(1316372988977975777/10411828384620436729956000000)*f2(n+2)*h+(105391490263288387/473264926573656214998000000)*f2(n+1)*h-(1284959669761615073/10411828384620436729956000000)*F2(n+2)*h^2+(72506125749079249/204153497737655622156000000)*h*f2(n))*u^12:
e4:=h^2*F2(n+1) = (60/23)*y2[n]-(60/23)*y2[n+1]+(25/46)*f2(n+2)*h+(32/23)*f2(n+1)*h-(4/23)*F2(n+2)*h^2+(31/46)*h*f2(n)+((209/3703)*y2[n]-(209/3703)*y2[n+1]+(1313/222180)*f2(n+2)*h+(1304/55545)*f2(n+1)*h-(131/18515)*F2(n+2)*h^2+(6011/222180)*h*f2(n))*u^2+((77491/35770980)*y2[n]-(77491/35770980)*y2[n+1]+(574843/2146258800)*f2(n+2)*h+(113536/134141175)*f2(n+1)*h-(53461/178854900)*F2(n+2)*h^2+(2258041/2146258800)*h*f2(n))*u^4+((151508243/1900512167400)*y2[n]-(151508243/1900512167400)*y2[n+1]+(1290306599/114030730044000)*f2(n+2)*h+(18919693/647901875250)*f2(n+1)*h-(113769323/9502560837000)*F2(n+2)*h^2+(4470322013/114030730044000)*h*f2(n))*u^6+((42120775181/14464625332248000)*y2[n]-(42120775181/14464625332248000)*y2[n+1]+(332746636891/734357901483360000)*f2(n+2)*h+(302396120633/298332897477615000)*f2(n+1)*h-(369019384141/795554393273640000)*F2(n+2)*h^2+(13797329479621/9546652719283680000)*h*f2(n))*u^8+((18953368786273/177347433208231440000)*y2[n]-(18953368786273/177347433208231440000)*y2[n+1]+(2430202319484337/138330997902420523200000)*f2(n+2)*h+(310803544671199/8645687368901282700000)*f2(n+1)*h-(203453960588449/11527583158535043600000)*F2(n+2)*h^2+(7380568619069419/138330997902420523200000)*h*f2(n))*u^10+((16436168060905785763/4164731353848174691982400000)*y2[n]-(16436168060905785763/4164731353848174691982400000)*y2[n+1]+(167160345356705269819/249883881230890481518944000000)*f2(n+2)*h+(461636091223370027/354948694930242161248500000)*f2(n+1)*h-(13852288092290788813/20823656769240873459912000000)*F2(n+2)*h^2+(29059878239787610409/14699051837111204795232000000)*h*f2(n))*u^12:
e5:=y3[n+2] = (7/23)*y3[n]+(16/23)*y3[n+1]+(12/23)*f3(n+2)*h+(16/23)*f3(n+1)*h-(2/23)*F3(n+2)*h^2+(2/23)*h*f3(n)+((24/3703)*y3[n]-(24/3703)*y3[n+1]+(48/18515)*f3(n+2)*h+(8/55545)*f3(n+1)*h-(116/55545)*F3(n+2)*h^2+(208/55545)*h*f3(n))*u^2+((901/2980915)*y3[n]-(901/2980915)*y3[n+1]+(7109/89427450)*f3(n+2)*h+(923/14904575)*f3(n+1)*h-(6241/89427450)*F3(n+2)*h^2+(14383/89427450)*h*f3(n))*u^4+((1979723/158376013950)*y3[n]-(1979723/158376013950)*y3[n+1]+(6364571/2375640209250)*f3(n+2)*h+(728327/215967291750)*f3(n+1)*h-(11785633/4751280418500)*F3(n+2)*h^2+(5106559/791880069750)*h*f3(n))*u^6+((6488435581/13259239887894000)*y3[n]-(6488435581/13259239887894000)*y3[n+1]+(8693517709/91794737685420000)*f3(n+2)*h+(260601208141/1789997384865690000)*f3(n+1)*h-(323357994149/3579994769731380000)*F3(n+2)*h^2+(891627999937/3579994769731380000)*h*f3(n))*u^8+((25090513463/1343541160668420000)*y3[n]-(25090513463/1343541160668420000)*y3[n+1]+(190450718149/55421072877572325000)*f3(n+2)*h+(47563947061/8210529315195900000)*f3(n+1)*h-(1475729910283/443368583020578600000)*F3(n+2)*h^2+(261738159769/27710536438786162500)*h*f3(n))*u^10+((244426606265778733/347060946154014557665200000)*y3[n]-(244426606265778733/347060946154014557665200000)*y3[n+1]+(1316372988977975777/10411828384620436729956000000)*f3(n+2)*h+(105391490263288387/473264926573656214998000000)*f3(n+1)*h-(1284959669761615073/10411828384620436729956000000)*F3(n+2)*h^2+(72506125749079249/204153497737655622156000000)*h*f3(n))*u^12:
e6:=h^2*F3(n+1) = (60/23)*y3[n]-(60/23)*y3[n+1]+(25/46)*f3(n+2)*h+(32/23)*f3(n+1)*h-(4/23)*F3(n+2)*h^2+(31/46)*h*f3(n)+((209/3703)*y3[n]-(209/3703)*y3[n+1]+(1313/222180)*f3(n+2)*h+(1304/55545)*f3(n+1)*h-(131/18515)*F3(n+2)*h^2+(6011/222180)*h*f3(n))*u^2+((77491/35770980)*y3[n]-(77491/35770980)*y3[n+1]+(574843/2146258800)*f3(n+2)*h+(113536/134141175)*f3(n+1)*h-(53461/178854900)*F3(n+2)*h^2+(2258041/2146258800)*h*f3(n))*u^4+((151508243/1900512167400)*y3[n]-(151508243/1900512167400)*y3[n+1]+(1290306599/114030730044000)*f3(n+2)*h+(18919693/647901875250)*f3(n+1)*h-(113769323/9502560837000)*F3(n+2)*h^2+(4470322013/114030730044000)*h*f3(n))*u^6+((42120775181/14464625332248000)*y3[n]-(42120775181/14464625332248000)*y3[n+1]+(332746636891/734357901483360000)*f3(n+2)*h+(302396120633/298332897477615000)*f3(n+1)*h-(369019384141/795554393273640000)*F3(n+2)*h^2+(13797329479621/9546652719283680000)*h*f3(n))*u^8+((18953368786273/177347433208231440000)*y3[n]-(18953368786273/177347433208231440000)*y3[n+1]+(2430202319484337/138330997902420523200000)*f3(n+2)*h+(310803544671199/8645687368901282700000)*f3(n+1)*h-(203453960588449/11527583158535043600000)*F3(n+2)*h^2+(7380568619069419/138330997902420523200000)*h*f3(n))*u^10+((16436168060905785763/4164731353848174691982400000)*y3[n]-(16436168060905785763/4164731353848174691982400000)*y3[n+1]+(167160345356705269819/249883881230890481518944000000)*f3(n+2)*h+(461636091223370027/354948694930242161248500000)*f3(n+1)*h-(13852288092290788813/20823656769240873459912000000)*F3(n+2)*h^2+(29059878239787610409/14699051837111204795232000000)*h*f3(n))*u^12:
e7:=y4[n+2] = (7/23)*y4[n]+(16/23)*y4[n+1]+(12/23)*f4(n+2)*h+(16/23)*f4(n+1)*h-(2/23)*F4(n+2)*h^2+(2/23)*h*f4(n)+((24/3703)*y4[n]-(24/3703)*y4[n+1]+(48/18515)*f4(n+2)*h+(8/55545)*f4(n+1)*h-(116/55545)*F4(n+2)*h^2+(208/55545)*h*f4(n))*u^2+((901/2980915)*y4[n]-(901/2980915)*y4[n+1]+(7109/89427450)*f4(n+2)*h+(923/14904575)*f4(n+1)*h-(6241/89427450)*F4(n+2)*h^2+(14383/89427450)*h*f4(n))*u^4+((1979723/158376013950)*y4[n]-(1979723/158376013950)*y4[n+1]+(6364571/2375640209250)*f4(n+2)*h+(728327/215967291750)*f4(n+1)*h-(11785633/4751280418500)*F4(n+2)*h^2+(5106559/791880069750)*h*f4(n))*u^6+((6488435581/13259239887894000)*y4[n]-(6488435581/13259239887894000)*y4[n+1]+(8693517709/91794737685420000)*f4(n+2)*h+(260601208141/1789997384865690000)*f4(n+1)*h-(323357994149/3579994769731380000)*F4(n+2)*h^2+(891627999937/3579994769731380000)*h*f4(n))*u^8+((25090513463/1343541160668420000)*y4[n]-(25090513463/1343541160668420000)*y4[n+1]+(190450718149/55421072877572325000)*f4(n+2)*h+(47563947061/8210529315195900000)*f4(n+1)*h-(1475729910283/443368583020578600000)*F4(n+2)*h^2+(261738159769/27710536438786162500)*h*f4(n))*u^10+((244426606265778733/347060946154014557665200000)*y4[n]-(244426606265778733/347060946154014557665200000)*y4[n+1]+(1316372988977975777/10411828384620436729956000000)*f4(n+2)*h+(105391490263288387/473264926573656214998000000)*f4(n+1)*h-(1284959669761615073/10411828384620436729956000000)*F4(n+2)*h^2+(72506125749079249/204153497737655622156000000)*h*f4(n))*u^12:
e8:=h^2*F4(n+1) = (60/23)*y4[n]-(60/23)*y4[n+1]+(25/46)*f4(n+2)*h+(32/23)*f4(n+1)*h-(4/23)*F4(n+2)*h^2+(31/46)*h*f4(n)+((209/3703)*y4[n]-(209/3703)*y4[n+1]+(1313/222180)*f4(n+2)*h+(1304/55545)*f4(n+1)*h-(131/18515)*F4(n+2)*h^2+(6011/222180)*h*f4(n))*u^2+((77491/35770980)*y4[n]-(77491/35770980)*y4[n+1]+(574843/2146258800)*f4(n+2)*h+(113536/134141175)*f4(n+1)*h-(53461/178854900)*F4(n+2)*h^2+(2258041/2146258800)*h*f4(n))*u^4+((151508243/1900512167400)*y4[n]-(151508243/1900512167400)*y4[n+1]+(1290306599/114030730044000)*f4(n+2)*h+(18919693/647901875250)*f4(n+1)*h-(113769323/9502560837000)*F4(n+2)*h^2+(4470322013/114030730044000)*h*f4(n))*u^6+((42120775181/14464625332248000)*y4[n]-(42120775181/14464625332248000)*y4[n+1]+(332746636891/734357901483360000)*f4(n+2)*h+(302396120633/298332897477615000)*f4(n+1)*h-(369019384141/795554393273640000)*F4(n+2)*h^2+(13797329479621/9546652719283680000)*h*f4(n))*u^8+((18953368786273/177347433208231440000)*y4[n]-(18953368786273/177347433208231440000)*y4[n+1]+(2430202319484337/138330997902420523200000)*f4(n+2)*h+(310803544671199/8645687368901282700000)*f4(n+1)*h-(203453960588449/11527583158535043600000)*F4(n+2)*h^2+(7380568619069419/138330997902420523200000)*h*f4(n))*u^10+((16436168060905785763/4164731353848174691982400000)*y4[n]-(16436168060905785763/4164731353848174691982400000)*y4[n+1]+(167160345356705269819/249883881230890481518944000000)*f4(n+2)*h+(461636091223370027/354948694930242161248500000)*f4(n+1)*h-(13852288092290788813/20823656769240873459912000000)*F4(n+2)*h^2+(29059878239787610409/14699051837111204795232000000)*h*f4(n))*u^12:
# Display of the solutions
h:=evalf(Pi/6):
omega:=1.0:
u:=omega*h:
N:=solve(h*p = 12*Pi/6, p):
n:=0:
exy1:= [seq](eval(cos(i)+0.0005*i*sin(i)), i=h..N,h):
exy2:= [seq](eval(-0.9995*sin(i)+0.0005), i=h..N,h):
exy3:= [seq](eval(sin(i)-0.0005*i*cos(i)), i=h..N,h):
exy4:= [seq](eval(0.9995*sin(i)+0.0005*i*sin(i)), i=h..N,h):
iny1:=1:
iny2:=0:
iny3:=0:
iny4:=0.9995:
err1 := Vector(N):
err2 := Vector(N):
c:=1:
inx:=0:
vars := y1[n+1],y1[n+2],y2[n+1],y2[n+2],y3[n+1],y3[n+2],y4[n+1],y4[n+2]:
for j from 0 to 2 do
x[j]:=inx+j*h:
end do:
printf("%4s%9s%9s%9s%9s%9s%9s%10s%10s%9s%9s%9s%10s\n",
"h","numy1","numy2","numy3","numy4",
"exy1","exy2","exy3","exy4",
"erry1","erry2","erry3","erry4");
st := time():
for k from 1 to N/2 do
param1:=y1[n]=iny1,y2[n]=iny2,y3[n]=iny3,y4[n]=iny4:
param2:=t[n]=x[0],t[n+1]=x[1],t[n+2]=x[2]:
res:=eval(<vars>, fsolve(eval({e||(1..8)},[param1,param2]),{vars})):
for i from 1 to 2 do
printf("%5.2f%9.3f%9.3f%9.3f%9.3f %8.5f%10.5f%10.5f%10.5f %8.2g%9.3g%9.3g%8.3g\n",
h*c,res[i],res[i+2],res[i+4],res[i+6],
exy1[c],exy2[c],exy3[c],exy4[c],
abs(res[i]-exy1[c]),abs(res[i+2]-exy2[c]),abs(res[i+4]-exy3[c]),abs(res[i+6]-exy4[c])):
err1[c] := abs(evalf(res[i]-exy1)):
err2[c] := abs(evalf(res[i+4]-exy3)):
c:=c+1:
end do:
iny1:=res[2]:
iny2:=res[4]:
iny3:=res[6]:
iny4:=res[8]:
inx:=x[2]:
for j from 0 to 2 do
x[j]:=inx+j*h:
end do:
end do:
v:=time() - st;
printf("Maximum error is %.13g", max(err1));
printf("Maximum error is %.13g", max(err2));