forcedosc.m

一个关于数学建模的资料,是纽约大学的程序和文档

% This program solves the equation for a forced oscilator

%  X_tt + X = f(t), X(0)=X0, X_t(0)=Y0

% through a predictor-corrector methodology. First, the
% equation is rewritten as a first order system:

%  X_t = Y
%  Y_t = -X + force(t).

% At time t, the present values of X(t) and Y(t), say X0 and Y0, 
% are used to compute
% a first estimate of the time derivatives X0_t and Y0_t. With
% these, a tentative value of X(t+dt) and Y(t+dt), say X1 and Y1,
% is computed through

%  X1 = X0 + dt*X0_t
%  Y1 = Y0 + dt*Y0_t.

% With X1, Y1 and t+dt, a new estimate for the time derivatives
% X1_t and Y1_t is computed. Finally, an approximate value of the
% solution at t+dt is given by

%  X(t+dt) = X(t) + dt/2 * (X0_t + X1_t)
%  Y(t+dt) = Y(t) + dt/2 * (Y0_t + Y1_t)

% t0, tf, dt, t, u (plotting variable):

t0=0; tf=150; dt=0.1; t=[t0:dt:tf]; u=zeros(size(t)); nt=size(t,2);

% Initial data:

X0=0; Y0=0;

% Period of the forcing:

P = 13.0;

u(1)=X0;

% Main loop over time:

for j=1:nt-1,

  X0t=Y0;
  Y0t=-X0+force(t(j),P);

  X1=X0+dt*X0t;
  Y1=Y0+dt*Y0t;

  X1t=Y1;
  Y1t=-X1+force(t(j+1),P);

  X0=X0+0.5*dt*(X0t+X1t);
  Y0=Y0+0.5*dt*(Y0t+Y1t);

  u(j+1)=X0;

end

% Plot results:

plot(t,u); xlabel('t'); ylabel('x');