restart;
with(plots): with(LinearAlgebra):
# TFSB Coefficients (symbolic in u)
beta0 := u -> (sin(u)*u^3 - 12*u^2 - 24*cos(u) + 24)/(12*(sin(u)*u + 2*cos(u) - 2)*u^2):
beta1 := u -> (5*sin(u)*u^3 + 12*cos(u)*u^2 + 24*cos(u) - 24)/(6*(sin(u)*u + 2*cos(u) - 2)*u^2):
beta2 := u -> beta0(u):
rho0 := u -> ((-u^2-12)*cos(u) - 5*u^2 + 12)/(12*(sin(u)*u + 2*cos(u) - 2)*u^2):
rho1 := u -> (-7*cos(u)*u^3 + 27*sin(u)*u^2 + 120*sin(u) - 120*u)/(60*u^2*(cos(u)*u + 2*u - 3*sin(u))):
rho2 := u -> -rho0(u):
# Secondary coefficients (simplified versions)
beta00 := u -> 13/42 - 9*u^2/7840:
beta10 := u -> 1/6 + u^2/720:
beta20 := u -> 1/42 - 17*u^2/70560:
beta01 := u -> 187/1680 + 611*u^2/705600:
beta11 := u -> 11/30 - 29*u^2/25200:
beta21 := u -> 37/1680 + 67*u^2/235200:
beta02 := u -> 11/70 + 491*u^2/352800:
beta12 := u -> 9/10 - 31*u^2/8400:
beta22 := u -> 31/70 + 811*u^2/352800:
rho01 := u -> 2/105 + 407*u^2/1058400:
rho11 := u -> -19/210 + 41*u^2/105840:
rho21 := u -> -1/168 - 101*u^2/529200:
rho02 := u -> 53/1680 + 1633*u^2/2116800:
rho12 := u -> 8/105 - 4*u^2/6615:
rho22 := u -> -101/1680 - 2273*u^2/2116800:
# Problem definition
omega := 1:
epsilon := 3*Pi/2:
phi := x -> 3*sin(x) - 5*cos(x): # history function
f := (x, v, vp, vd) -> -v - vd + 3*cos(x) + 5*sin(x):
g := proc(x, v, vp, vd, vdp)
local fx, fv, fvp, fvd;
fx := -3*sin(x) + 5*cos(x);
fv := -1;
fvp := 0;
fvd := -1;
return fx + fv*vp + fvp*0 + fvd*vdp;
end proc:
# Initial conditions
a := 0: b := 10:
v0 := -5: vp0 := 3:
# Variable step-size parameters
tol := 1e-10:
h_min := 0.01:
h_max := 0.5:
h_init := Pi/8:
# Store results
X := [a]: V := [v0]: Vp := [vp0]:
h_curr := h_init:
x_curr := a:
v_curr := v0:
vp_curr := vp0:
# For history: need v at x-epsilon
get_v_delayed := proc(xx)
if xx < a then return phi(xx);
else
# Interpolate from stored solution
idx := 1;
while idx < nops(X) and X[idx] < xx do idx := idx+1; end do;
if idx = 1 then return phi(xx);
elif X[idx] = xx then return V[idx];
else
# Linear interpolation
return V[idx-1] + (V[idx]-V[idx-1])*(xx-X[idx-1])/(X[idx]-X[idx-1]);
end if;
end if;
end proc:
# Newton solver for block
solve_block := proc(x0, v0, vp0, h, omega)
local u, bet0, bet1, bet2, rho0, rho1, rho2, bet00, bet10, bet20, bet01, bet11, bet21, bet02, bet12, bet22,
rho01, rho11, rho21, rho02, rho12, rho22, F, J, V0, V1, V2, Vp0, Vp1, Vp2, tolN, iter, dv, dV;
u := omega*h;
bet0 := beta0(u); bet1 := beta1(u); bet2 := beta2(u);
rho0 := rho0(u); rho1 := rho1(u); rho2 := rho2(u);
bet00 := beta00(u); bet10 := beta10(u); bet20 := beta20(u);
bet01 := beta01(u); bet11 := beta11(u); bet21 := beta21(u);
bet02 := beta02(u); bet12 := beta12(u); bet22 := beta22(u);
rho01 := rho01(u); rho11 := rho11(u); rho21 := rho21(u);
rho02 := rho02(u); rho12 := rho12(u); rho22 := rho22(u);
# Initial guesses
V1 := v0 + h*vp0;
V2 := v0 + 2*h*vp0;
Vp1 := vp0;
Vp2 := vp0;
tolN := 1e-12;
for iter from 1 to 10 do
# Compute delayed values
vd0 := get_v_delayed(x0 - epsilon);
vd1 := get_v_delayed(x0 + h - epsilon);
vd2 := get_v_delayed(x0 + 2*h - epsilon);
vdp0 := (get_v_delayed(x0 - epsilon + 1e-8) - vd0)/1e-8;
vdp1 := (get_v_delayed(x0 + h - epsilon + 1e-8) - vd1)/1e-8;
vdp2 := (get_v_delayed(x0 + 2*h - epsilon + 1e-8) - vd2)/1e-8;
# Compute gamma and g
gam0 := f(x0, v0, vp0, vd0);
gam1 := f(x0+h, V1, Vp1, vd1);
gam2 := f(x0+2*h, V2, Vp2, vd2);
g0 := g(x0, v0, vp0, vd0, vdp0);
g1 := g(x0+h, V1, Vp1, vd1, vdp1);
g2 := g(x0+2*h, V2, Vp2, vd2, vdp2);
# Residuals
F1 := h*vp0 - (V1 - v0 + h^2*(bet00*gam0 + bet10*gam1 + bet20*gam2)
+ h^3*(rho01*g0 + rho11*g1 + rho21*g2));
F2 := h*Vp1 - (V1 - v0 + h^2*(bet01*gam0 + bet11*gam1 + bet21*gam2)
+ h^3*(rho01*g0 + rho11*g1 + rho21*g2));
F3 := h*Vp2 - (V1 - v0 + h^2*(bet02*gam0 + bet12*gam1 + bet22*gam2)
+ h^3*(rho02*g0 + rho12*g1 + rho22*g2));
F4 := V2 - (2*V1 - v0 + h^2*(bet0*gam0 + bet1*gam1 + bet2*gam2)
+ h^3*(rho0*g0 + rho1*g1 + rho2*g2));
F := Vector([F1, F2, F3, F4]);
if LinearAlgebra:-Norm(F) < tolN then break; end if;
# Approximate Jacobian (finite differences)
J := Matrix(4,4);
delta := 1e-6;
for j from 1 to 4 do
V_pert := Vector([V1, V2, Vp1, Vp2]);
V_pert[j] := V_pert[j] + delta;
V1p := V_pert[1]; V2p := V_pert[2]; Vp1p := V_pert[3]; Vp2p := V_pert[4];
gam1p := f(x0+h, V1p, Vp1p, get_v_delayed(x0+h-epsilon));
gam2p := f(x0+2*h, V2p, Vp2p, get_v_delayed(x0+2*h-epsilon));
g1p := g(x0+h, V1p, Vp1p, get_v_delayed(x0+h-epsilon),
(get_v_delayed(x0+h-epsilon+1e-8)-get_v_delayed(x0+h-epsilon))/1e-8);
g2p := g(x0+2*h, V2p, Vp2p, get_v_delayed(x0+2*h-epsilon),
(get_v_delayed(x0+2*h-epsilon+1e-8)-get_v_delayed(x0+2*h-epsilon))/1e-8);
F1p := h*vp0 - (V1p - v0 + h^2*(bet00*gam0 + bet10*gam1p + bet20*gam2p)
+ h^3*(rho01*g0 + rho11*g1p + rho21*g2p));
F2p := h*Vp1p - (V1p - v0 + h^2*(bet01*gam0 + bet11*gam1p + bet21*gam2p)
+ h^3*(rho01*g0 + rho11*g1p + rho21*g2p));
F3p := h*Vp2p - (V1p - v0 + h^2*(bet02*gam0 + bet12*gam1p + bet22*gam2p)
+ h^3*(rho02*g0 + rho12*g1p + rho22*g2p));
F4p := V2p - (2*V1p - v0 + h^2*(bet0*gam0 + bet1*gam1p + bet2*gam2p)
+ h^3*(rho0*g0 + rho1*g1p + rho2*g2p));
Fp := Vector([F1p, F2p, F3p, F4p]);
J[1..4, j] := (Fp - F)/delta;
end do;
dV := LinearAlgebra:-LinearSolve(J, -F);
V1 := V1 + dV[1]; V2 := V2 + dV[2]; Vp1 := Vp1 + dV[3]; Vp2 := Vp2 + dV[4];
end do;
return [V1, V2, Vp1, Vp2];
end proc:
# Main variable step-size loop
printf("Variable step-size integration for Example 1\n");
printf("tol = %e, h_init = %f\n", tol, h_init);
while x_curr < b - 1e-12 do
# Try current step
sol := solve_block(x_curr, v_curr, vp_curr, h_curr, omega);
V1 := sol[1]; V2 := sol[2]; Vp1 := sol[3]; Vp2 := sol[4];
# Compute with two half-steps
sol_half1 := solve_block(x_curr, v_curr, vp_curr, h_curr/2, omega);
V_mid := sol_half1[2]; Vp_mid := sol_half1[4];
sol_half2 := solve_block(x_curr + h_curr/2, V_mid, Vp_mid, h_curr/2, omega);
V2_half := sol_half2[2];
# Error estimate
err := abs(V2 - V2_half) / (2^6 - 1);
if err < tol then
# Accept step
x_next := x_curr + 2*h_curr;
X := [op(X), x_curr + h_curr, x_next];
V := [op(V), V1, V2];
Vp := [op(Vp), Vp1, Vp2];
x_curr := x_next;
v_curr := V2;
vp_curr := Vp2;
printf("x = %7.4f, h = %8.5f, err = %12.5e\n", x_curr, h_curr, err);
# Adjust step size
if err < tol/2 then
h_curr := min(2*h_curr, h_max);
end if;
else
# Reject step, reduce h
h_curr := max(h_curr/2, h_min);
printf(" Rejecting, new h = %8.5f\n", h_curr);
end if;
end do:
# Exact solution for comparison
exact := x -> 3*sin(x) - 5*cos(x);
errors := [seq(abs(V[i] - exact(X[i])), i=1..nops(X))];
# Visualization
p1 := pointplot([seq([X[i], errors[i]], i=1..nops(X))], color=red, symbol=circle,
title="Example 1: Variable Step-Size TFSB - Absolute Errors",
labels=["x", "Error"], labeldirections=[horizontal,vertical]);
p2 := plot([[x_curr, h_curr]], x=a..b, style=point, color=blue,
title="Step-size evolution", labels=["x", "h"]);
display(p1);
display(p2);
printf("\nFinal results for Example 1:\n");
printf("Number of steps: %d\n", nops(X)-1);
printf("Maximum error: %e\n", max(errors));
printf("Final step-size: %f\n", h_curr);
