adamsbashforth.m

↓学部四年生のときにとった海洋モデリングの宿題。

%clear all;

dt=1;

dx=0.1;

dy=0.1;

du=0.1;

dv=0.1;

period=1000;

nmax=period/dt;

u0=0.05;

v0=0.05;

f=0.5*2*pi/24; % coriolis parameter at mid-latitude (i.e. at 45 deg.N)

u(1)=0.7;

v(1)=0;

x(1)=0;

y(1)=0;

%solve for u(2) and v(2) using forward scheme since we don't have it.

u(2)=u(1)*dt;

v(2)=v(1)*dt;

for n=2:nmax-1

u(n+1)=u(n)+f*dt*0.5*(3*v(n)-v(n-1));

v(n+1)=v(n)-f*dt*0.5*(3*u(n)-u(n-1));

end

for n=1:nmax-1

%euler forward

x(n+1)=x(n)+dt*(u0+u(n));

y(n+1)=y(n)+dt*(v0+v(n));

end

figure;plot(x,y);title 'Trajectory of the parcel';xlabel 'x(m)'; ylabel 'y(m)';