In which our hero asks for help.

Jun 18, 2008 15:40

I'm having trouble with a program I'm writing for my research project. It's telling me a I have things when I shouldn't. Read through the code and figure out what's wrong if you can. The problematic text is bolded, the not so problematic text is struck out.

function [kink]=crop(n,t,dt)
% Created June 18, 2008. A program to count the kinks of a annular system
% of coupled harmonic oscillators. A kink is defined rather difficultly.
% Mathematically:
% If a(i-1)-a(i) and a(i+1)-a(i) share the same sign, then there exists a
% kink at a(i)
% Where a(i) is the phase angle of oscillator i. There are special
% situations where a(i±1)=a(i) or a(i±1)=a(i)±? which need to be discussed
% in some detail (080618)
%
% For random inital conditions
%
% Inputs
%
% n is the number of oscillators
%
% t is the duration of the event
%
% dt is the step size between t's

% m is the number of t terms:
m = t/dt;

%% Initialising the Arrays and Initial Conditions

a = zeros(n,m+1);

% The terms here are randomly generated
a(:,1) = 2*pi*rand(1,n);

%% Determining Initial Number of Kinks
% Added after harness.m

k = 0;

for r=1

a0 = a(n,1);
a1 = a(1,1);
a2 = a(2,1);

if sign(sin(a0-a1))~=0
if sign(sin(a0-a1))==sign(sin(a2-a1))

k=k+1;

end
end

end

for r=2:n-1

a0 = a(r-1,1);
a1 = a(r,1);
a2 = a(r+1,1);

if sign(sin(a0-a1))~=0
if sign(sin(a0-a1))==sign(sin(a2-a1))

k=k+1;

end
end

end

for r=n-1

a0 = a(n-1,1);
a1 = a(n,1);
a2 = a(1,1);

if sign(sin(a0-a1))~=0
if sign(sin(a0-a1))==sign(sin(a2-a1))

k=k+1;

end
end

end

kink(1)=k;

%% Generation of future terms
% Terms for b0, b1,b2 and k and kink were added after harness.m

for s=2:m+1

k = 0;

for r=1

a0 = a(n,s-1);
a1 = a(1,s-1);
a2 = a(2,s-1);

da = sin(a0-a1) + sin(a2-a1);

a(r,s) = mod(a1+da*dt,2*pi);

b0 = a(n,s);
b1 = a(1,s);
b2 = a(2,s);

if sign(sin(b0-b1))~=0
if sign(sin(b0-b1))==sign(sin(b2-b1))

k=k+1;

end
end

end

for r=2:(n-1)

a0 = a(r-1,s-1);
a1 = a(r,s-1);
a2 = a(r+1,s-1);

da = sin(a0-a1) + sin(a2-a1);

a(r,s) = mod(a1+da*dt,2*pi);

b0 = a(r-1,s);
b1 = a(r,s);
b2 = a(r+1,s);

if sign(sin(b0-b1))~=0
if sign(sin(b0-b1))==sign(sin(b2-b1))

k=k+1;

end
end

end

for r=n

a0 = a(n-1,s-1);
a1 = a(n,s-1);
a2 = a(1,s-1);

da = sin(a0-a1) + sin(a2-a1);

a(r,s) = mod(a1+da*dt,2*pi);

b0 = a(n-1,s);
b1 = a(n,s);
b2 = a(1,s);

if sign(sin(b0-b1))~=0
if sign(sin(b0-b1))==sign(sin(b2-b1))

k=k+1;

end
end

end

kink(s) = k;

end

aend = zeros(n,2);
aend(:,1) = a(:,1);
aend(:,2) = a(:,m+1);
aend

%% Graphically portraying the solution

N = 1:n;
T = 0:dt:t;

figure % Added after harness.m
mesh(T,N,a)

%% Plotting the Number of Kinks as a Function of Time
% added after harness.m

figure
plot(T,kink)

%% End Program
Previous post Next post
Up