Svår matlab uppgift. Behöver hjälp att komma igång!

sitter med ett svårt matlab problem och har ingen ide om hur man ska börja. är det någon som kan ge mig några pointers på vad man kan göra.

det jag tror mig veta är att man skall använda ode45.

Läsa detta kan vara en bra början.

http://www.mathworks.com/help/techdoc/ref/ode23.html

Är det något specifikt du inte förstår så specificera gärna vad det är.

Här kan du få ett exempel på hur filen där du definierar differentialekvationen ska se ut. Detta exempel är från någon gång för länge sedan då jag roade mig med att räkna på en railgun med koaxialt tvärsnitt.

Funktionen som sparas som diffsystem.m:
function xprime = diffsystem(t,x,PAR)

%test for the solution
% b=20e-3;
% a=10e-3;
% L=1e-7;
mu0=4*pi*1e-7;
% m=1e-3;
% Lprim=1e-7;
% C=1e-3;
% R=0.0001;

%PAR=[a b R C Lrail LprimRail]
a=PAR(1);
b=PAR(2);
R=PAR(3);
C=PAR(4);
Lrail=PAR(5);
LprimRail=PAR(6);
m=PAR(7);
L0=PAR(8)

%L=l*(mu0/(2*pi))*(0.25 + log(b/a)); %Total inductance of the rails

%LprimRail=Lrail/l; %Inductance gradient

%zrail=[0:l/100:l];
%Lz=LprimRail*zrail;


%Lrail=interp1(zrail,Lz,x(length(x),4),'spline');

% Since the states are passed in as a single vector, let
% x(1) = Q =dq/dt = I
% x(2) = v =speed
% x(3) = q =charge
% x(4) = z =position


%Eq1
xprime(1)=-(1/(L0 + LprimRail*x(4)))*((1/C)*x(3) + x(1)*LprimRail*x(2) + R*x(1));
xprime(2)=mu0*((x(1))^2)*log(b/a)/(2*pi*m);
xprime(3)=x(1);
xprime(4)=x(2);


xprime = xprime(;
% This ensures that the vector returned is a column vector




Anropet sker från main_railgun.m:


InitCond=[Q0 v0 q0 z0]

PAR=[a b R C Lrail LprimRail m L0]

t0 = 0; % Start time
tf = 1e-2; % Stop time
tode=[0:tf/1000:tf]
options=[];
[t,x] = ode45(@diffsystem,tode,InitCond,options,PAR);
Q = x(:,1);
v = x(:,2);
q = x(:,3);
z = x(:,4);


Notera att jag har plockat bort alla deklarareringar/beräkningar av parametrar. Detta är enbart för att visa strukturen som ditt problem kommer att få.
Lycka till

edit. Det blev visst en smiley mitt i koden, men du förstår att det ska vara ( : )

Tja sitter med samma uppgift har löst första delen med att bara plotta partiklarnas flöde. Nu behöver jag hjälp med "Den effektiva algoritmen"

Såhär ser programmet ut för tillfället:
%Plottar partiklarnas läge
R=2;
x=[-4 -4 -4 -4];
y=[0.2 0.6 1.0 1.4];
endp=[];
for i=1:4;
ustart=[x(i) y(i)];
[t,U]=ode45('F',[0 12],ustart);
X=U(:,1)
Y=U(:,2)
plot(X,Y)
hold on
endp=[endp [X(end) Y(end)]']
end
Xny=[endp(2,1) endp(2,4)]
Yny=[endp(1,1) endp(1,4)]
for i=1:2
ustart=[Xny(i) Yny(i)];
[t,U]=ode45('F',[0 12],ustart);
X=U(:,1)
Y=U(:,2)
plot(X,Y)
end

%Plottar cylindern
theta=[0: pi/100: pi]; %För hel cirkel skriv 2*pi
X=R*cos(theta);
Y=R*sin(theta);
axis equal
plot(X,Y)

Sista for loopen kan ni bortse ifrån testade bara olika sätt att tackla problemet.

Har även en funktionsfil som ser ut på detta vis:

function prim=F(t,u)
R=2
dx=1-(R^2*(u(1)^2-u(2)^2))/((u(1)^2+u(2)^2)^2);
dy=-(2*u(1)*u(2)*(R^2))/(u(1)^2+u(2)^2)^2;
prim=[dx;dy]

Har gjort lite mellanrum på vissa ställen så det inte ska bli massa smileys