| > |
restart;
with(plots):
# ============================================================
# NONLINEAR WAVE EQUATION SOLVED AS COUPLED ODEs
# ============================================================
#
# ORIGINAL PDE: ∂²u/∂t² = c²(λ(H + ∂u/∂x)) · ∂²u/∂x²
#
# AFTER SPATIAL DISCRETIZATION (FINITE DIFFERENCES):
# We obtain a SYSTEM of COUPLED ORDINARY DIFFERENTIAL EQUATIONS (ODEs):
#
# For each internal point i = 1, 2, ..., Nx-1:
#
# d²u_i/dt² = c²(λ(H + (u_{i+1}-u_{i-1})/(2Δx))) · (u_{i+1}-2u_i+u_{i-1})/Δx²
#
# This is a system of (Nx-1) coupled, nonlinear ODEs.
# ============================================================
# ============================================================
# STEP 1: PHYSICAL SYSTEM PARAMETERS
# ============================================================
L := 0.04; # Length of the rod [m] (4 cm)
A0 := 12*L/2*10^(-6); # Initial amplitude [m] (~2.4e-7 m)
Nx := 100; # Number of spatial points (creates 99 ODEs)
dx := L/Nx; # Spatial step size [m]
x_vals := [seq(i*dx, i=0..Nx)]; # Positions along the rod [m]
printf("========================================\n");
printf(" COUPLED ODEs FROM WAVE EQUATION\n");
printf("========================================\n\n");
printf("STEP 1: SYSTEM PARAMETERS\n");
printf(" Rod length L = %.4f m (%.1f cm)\n", L, L*100);
printf(" Initial amplitude A0 = %.2e m\n", A0);
printf(" Number of spatial points Nx+1 = %d\n", Nx+1);
printf(" Number of COUPLED ODEs = %d (for internal points)\n", Nx-1);
printf(" Spatial step dx = %.2e m\n\n", dx);
# ============================================================
# STEP 2: MATERIAL MODEL (PIECEWISE FUNCTIONS)
# ============================================================
# c_sq(H) = squared speed of sound [m²/s²]
c_sq := proc(H)
if H < 700 then
return -3143*H + 1.99e7; # Region 1: decreasing speed
elif H < 1000 then
return 7333*H + 1.257e7; # Region 2: increasing speed
else
return 1.99e7; # Region 3: constant speed
end if;
end proc:
# lambda(H) = helper function for inverse relation
lambda := proc(H)
if H < 1000 then
return (3/250000000)*H;
else
return 3/250000;
end if;
end proc:
printf("STEP 2: MATERIAL MODEL\n");
printf(" c²(H) = piecewise function:\n");
printf(" H < 700: c² = -3143·H + 1.99e7\n");
printf(" 700 ≤ H < 1000: c² = 7333·H + 1.257e7\n");
printf(" H ≥ 1000: c² = 1.99e7 (constant)\n");
printf(" λ(H) = 3/250000000 * H for H < 1000\n\n");
# ============================================================
# STEP 3: NUMERICAL PARAMETERS
# ============================================================
c_max := sqrt(1.99e7); # Maximum wave speed [m/s]
dt := 0.5 * dx / c_max; # Time step (CFL = 0.5 for stability)
f0 := sqrt(c_sq(650))/(2*L); # Resonance frequency [Hz]
T0 := 1/f0; # Period [s]
Nt := 1500; # Number of time steps
printf("STEP 3: NUMERICAL PARAMETERS\n");
printf(" Maximum wave speed c_max = %.1f m/s\n", c_max);
printf(" Resonance frequency f0 = %.1f Hz (ultrasonic!)\n", f0);
printf(" Period T0 = %.2e s\n", T0);
printf(" Time step dt = %.2e s\n", dt);
printf(" Number of time steps Nt = %d\n", Nt);
printf(" CFL number = %.3f (must be < 1 for stability)\n\n", c_max*dt/dx);
# ============================================================
# STEP 4: INITIALIZE THE SOLUTION (Three time layers)
# ============================================================
# We need three arrays for the time integration:
# u_old = solution at time t - dt
# u_curr = solution at time t
# u_new = solution at time t + dt
u_old := Array(0..Nx);
u_curr := Array(0..Nx);
u_new := Array(0..Nx);
# INITIAL CONDITION: u(x,0) = A0·cos(πx/L)
# This is the FIRST HARMONIC vibration mode
for i from 0 to Nx do
u_curr[i] := A0 * cos(Pi * x_vals[i+1] / L);
u_old[i] := u_curr[i]; # Initial velocity = 0 => u_old = u_curr
end do:
printf("STEP 4: INITIAL CONDITION\n");
printf(" u(x,0) = A0·cos(πx/L)\n");
printf(" ∂u/∂t(x,0) = 0\n");
printf(" This creates the first harmonic vibration mode.\n");
printf(" The rod is maximally stretched at the middle (x = %.2f m).\n\n", L/2);
# Display the initial values for the first few points (educational)
printf(" Initial values for the first 5 points (i=0 to 4):\n");
for i from 0 to 4 do
printf(" u[%d](0) = %.2e m at x = %.4f m\n", i, u_curr[i], x_vals[i+1]);
end do;
printf(" ...\n\n");
# ============================================================
# STEP 5: BOUNDARY CONDITIONS (Neumann - Free Ends)
# ============================================================
# Neumann: ∂u/∂x = 0 at both ends
# This means: u[0] = u[1] and u[Nx] = u[Nx-1]
# Physically: The ends are FREE (no force)
apply_bc := proc()
u_new[0] := u_new[1];
u_new[Nx] := u_new[Nx-1];
end proc:
printf("STEP 5: BOUNDARY CONDITIONS\n");
printf(" Neumann: ∂u/∂x = 0 at x = 0 and x = L\n");
printf(" This means: u₀(t) = u₁(t) and u_%d(t) = u_%d(t)\n", Nx, Nx-1);
printf(" Physically: The ends are FREE (can move without force).\n");
printf(" Waves reflect at the ends WITHOUT sign reversal.\n\n");
# ============================================================
# STEP 6: THE COUPLED ODE SYSTEM - EXPLANATION
# ============================================================
printf("STEP 6: THE COUPLED ODE SYSTEM\n");
printf(" After spatial discretization, we get %d coupled ODEs:\n\n", Nx-1);
for i from 1 to 3 do
printf(" ODE %d: d²u_%d/dt² = c²(λ(H + (u_%d-u_%d)/(2dx))) · (u_%d-2u_%d+u_%d)/dx²\n",
i, i, i+1, i-1, i+1, i, i-1);
end do;
printf(" ...\n");
for i from Nx-3 to Nx-1 do
printf(" ODE %d: d²u_%d/dt² = c²(λ(H + (u_%d-u_%d)/(2dx))) · (u_%d-2u_%d+u_%d)/dx²\n",
i, i, i+1, i-1, i+1, i, i-1);
end do;
printf("\n Each ODE connects point i with its neighbors i-1 and i+1.\n");
printf(" This is why they are called COUPLED ODEs.\n\n");
# ============================================================
# STEP 7: SOLVING THE COUPLED ODEs USING FINITE DIFFERENCES
# ============================================================
# We use an explicit time integration scheme:
# u_new[i] = 2*u_curr[i] - u_old[i] + dt² · (d²u_i/dt²)
#
# This comes from the central difference approximation:
# d²u_i/dt² ≈ (u_new[i] - 2*u_curr[i] + u_old[i]) / dt²
printf("STEP 7: SOLVING THE COUPLED ODEs\n");
printf(" Using explicit time integration (central differences):\n");
printf(" u_new[i] = 2*u_curr[i] - u_old[i] + dt² · (d²u_i/dt²)\n\n");
printf(" Starting numerical integration...\n\n");
# Storage indices for visualization
save_indices := [42, 106, 211, 317, 422];
solutions := table();
# TIME INTEGRATION LOOP
for n from 1 to Nt do
# --------------------------------------------------------
# COMPUTE ALL COUPLED ODEs
# Each iteration of this loop computes ONE ODE for point i
# The ODEs are COUPLED because u_new[i] depends on u_curr[i-1] and u_curr[i+1]
# --------------------------------------------------------
for i from 1 to Nx-1 do
# ----- COUPLING TERM 1: LOCAL STRAIN (∂u/∂x) -----
# This connects u_i with its neighbors u_{i-1} and u_{i+1}
du_dx := (u_curr[i+1] - u_curr[i-1]) / (2 * dx);
# ----- MATERIAL PARAMETER H DEPENDS ON STRAIN -----
H_total := 650 + du_dx;
# ----- SQUARED SPEED OF SOUND c² -----
if H_total < 700 then
c2 := -3143*H_total + 1.99e7;
elif H_total < 1000 then
c2 := 7333*H_total + 1.257e7;
else
c2 := 1.99e7;
end if;
# ----- COUPLING TERM 2: SECOND DERIVATIVE (∂²u/∂x²) -----
# This also connects u_i with its neighbors u_{i-1} and u_{i+1}
u_xx := (u_curr[i+1] - 2*u_curr[i] + u_curr[i-1]) / dx^2;
# ----- THE COUPLED ODE FOR POINT i -----
# d²u_i/dt² = c² · u_xx
# Discretized in time:
u_new[i] := 2*u_curr[i] - u_old[i] + dt^2 * c2 * u_xx;
end do;
# --------------------------------------------------------
# END OF COUPLED ODE COMPUTATION
# --------------------------------------------------------
# Apply boundary conditions (coupling the end points)
apply_bc();
# Store solutions for visualization
for idx from 1 to 5 do
if n = save_indices[idx] then
solutions[n] := copy(u_new);
printf(" Stored: t = %.3f T0 (time step %d/%d)\n", n*dt/T0, n, Nt);
end if;
end do;
# Shift time layers for next iteration
for i from 0 to Nx do
u_old[i] := u_curr[i];
u_curr[i] := u_new[i];
end do;
# Progress indicator
if n mod 500 = 0 then
printf(" Progress: %d/%d (%.0f%%)\n", n, Nt, 100*n/Nt);
end if;
end do:
printf("\nSimulation completed!\n\n");
# ============================================================
# STEP 8: VISUALIZATION OF THE COUPLED ODE SOLUTION
# ============================================================
printf("STEP 8: VISUALIZATION\n");
printf(" Generating plots of the solution u_i(t)...\n");
colors := ["red", "blue", "green", "magenta", "black"];
plot_list := [];
for idx from 1 to 5 do
n_idx := save_indices[idx];
t_val := evalf(n_idx * dt / T0, 3);
points := [seq([x_vals[i+1], solutions[n_idx][i]], i=0..Nx)];
p := plot(points,
labels=["x [m]", "u(x,t) [m]"],
title=sprintf("Wave at t = %.3f T₀", t_val),
color=colors[idx],
legend=sprintf("t = %.3f T₀", t_val),
thickness=2);
plot_list := [op(plot_list), p];
end do:
display(plot_list, title="Solution of the Coupled ODE System", size=[800,500]);
printf(" Plot shows the wave at 5 different times.\n\n");
# ============================================================
# STEP 9: ANIMATION OF THE COUPLED ODE SOLUTION
# ============================================================
printf("STEP 9: ANIMATION\n");
printf(" Generating animation of the moving wave...\n");
animation_frames := [];
t_labels := [0.1, 0.25, 0.5, 0.75, 1.0];
explanation := [
"t = 0.10 T₀: Wave begins moving right.\nThe wave deforms because c² depends on strain.",
"t = 0.25 T₀: Wave approaches the right end.\nAmplitude changes due to nonlinearity.",
"t = 0.50 T₀: Wave reaches the end and reflects.\n∂u/∂x=0 at the end (free boundary).",
"t = 0.75 T₀: Wave moves back left.\nWave shape is mirrored.",
"t = 1.00 T₀: Wave returns to initial condition.\nOne complete period completed."
];
for idx from 1 to 5 do
n_idx := save_indices[idx];
t_val := evalf(n_idx * dt / T0, 3);
points := [seq([x_vals[i+1], solutions[n_idx][i]], i=0..Nx)];
p := plot(points,
labels=["x [m]", "u(x,t) [m]"],
title=sprintf("Longitudinal wave at t = %.2f T₀", t_labels[idx]),
color="blue",
thickness=3,
view=[0..L, -A0*2..A0*2]);
text := plots:-textplot([0.02, A0*1.5, explanation[idx]],
font=["Arial", 10],
color="darkred",
align="LEFT");
animation_frames := [op(animation_frames), display(p, text)];
end do:
# Repeat frames for smooth animation
smooth_frames := [];
for rep from 1 to 4 do
for idx from 1 to 5 do
smooth_frames := [op(smooth_frames), animation_frames[idx]];
end do;
end do:
display(smooth_frames, insequence=true, size=[800,600]);
printf(" Animation started! (%d frames)\n\n", nops(smooth_frames));
# ============================================================
# STEP 10: SUMMARY FOR THE STUDENT
# ============================================================
printf("========================================\n");
printf(" SUMMARY\n");
printf("========================================\n\n");
printf("WHAT IS A COUPLED ODE SYSTEM?\n");
printf(" • We have %d unknown functions u_i(t), i = 1..%d\n", Nx-1, Nx-1);
printf(" • Each ODE connects u_i with its neighbors u_{i-1} and u_{i+1}\n");
printf(" • This is why they are called COUPLED ODEs\n\n");
printf("WHERE ARE THE ODEs IN THE CODE?\n");
printf(" • Inside the main time loop (for n from 1 to Nt)\n");
printf(" • In the inner loop (for i from 1 to Nx-1)\n");
printf(" • Each iteration computes ONE ODE for point i\n\n");
printf("WHAT DOES EACH ODE DO?\n");
printf(" • Computes d²u_i/dt² = c² · (u_{i+1}-2u_i+u_{i-1})/dx²\n");
printf(" • Updates u_new[i] = 2*u_curr[i] - u_old[i] + dt² · (d²u_i/dt²)\n\n");
printf("ANSWER TO THE ORIGINAL QUESTION:\n");
printf(" Q: Can pdsolve handle c²(λ(H + ∂u/∂x))?\n");
printf(" A: NO, pdsolve CANNOT handle this nonlinear PDE.\n");
printf(" SOLUTION: Convert to COUPLED ODEs and use finite differences.\n");
printf("========================================\n");
|
