Solving system of ODEs in Matlab

Question

I am trying to code an ODE system for $n=[-2,-1,0,1]$

$\frac{dx_n}{dt}=V_nx_n-x_{n+1}-x_{n-1}$

with some IC, say $[0 ~ 0~ 0~0]$: I have also a B.C at the two ends $n=-2$ and $n=1$, which I don't know how to implement at the moment, so to start I just tried to code with only I.Cs but having some problem. The code is

tic;
x0 = [0 0 0 0];
tspan = [0, 2];
V=-2.5;

eoms = @(t,x) [V*x(1)-x(2); V*x(2)-x(3)-x(1); V*x(3)-x(4)-x(2); V*x(4)-x(5)-x(3)];
[t, x] = ode45(eoms, tspan, x0)
toc
plot(t,x)

where x(1) corresponds to $x_{-2}$ , x(2) corresponds to $x_{-1}$ and so on, for the four equations corresponding to $n=-2,-1,0,1$. The term $x_{-3}$ has been ignored. The code gives some errors "Index exceeds matrix dimensions."

I have actually a more complicated R.H.S, and the description is a simplified one. I have slightly corrected the Eq.1 which was not correctly posted before. Do you think if I can correctly create a vector form of the R.H.S (with B.Cs embedded) then it is possible to use ode45? — AtoZ, Feb 20 '19 at 05:06
The term corresponding to x(5) is then a boundary term for this. — AtoZ, Feb 20 '19 at 05:09
the point is that when you call odeXX, the dimensions of the state vector and the vector of initial conditions must coincide. Your initial conditions have dimension 4. So, matlab tries to compute your function eoms for a 4-dim vector. However, you try to access the 5th element, which generates the error. In your case, you should pass the boundary value to the function as a parameter — Dmitry, Feb 20 '19 at 14:21
https://www.mathworks.com/help/optim/ug/passing-extra-parameters.html — Dmitry, Feb 20 '19 at 14:25
@Dmitry Thank you. I have actually tried to do it via RK4, by modifying a code: The code is: — AtoZ, Feb 21 '19 at 11:36
@Dmitry: The code works fine for the $x_{n+1}+x_{n-1}$ part: but when I try to add the term V (which is coded as): z(1+(numel(p)-2)/2:2+(numel(p)-2)/2)=[v*y(1+ceil(end/2)),v*y(2+ceil(end/2))] with diff(y,2)+2*y(2:end-1) I get problems. — AtoZ, Feb 21 '19 at 11:44
So effectively the code is only for $dx_n/dt=x_{n+1}+x_{n-1}$. I was not successful in adding the $V_n$ contribution, as I described above. Also to add that $V_n$ is only nonzero at sites 1,2. This is why I wrote z(1+(numel(p)-2)/2:2+(numel(p)-2)/2)=[v*y(1+ceil(end/2)),v*y(2+ceil(end/2))] — AtoZ, Feb 21 '19 at 11:47
The data psi1,psi2,psi3,psi0,R0,R1 are just related to definitions: y is the actual (trial) solution which will be plugged into the Eq.1 (i.e., the $x$s). I wrote the vector for V_n (z(1+(numel(p)-2)/2:2+(numel(p)-2)/2)=[v*y(1+ceil(end/2)),v*y(2+ceil(end/2))])correctly, but could not add it to diff(y,2)+2*y(2:end-1) ($x_{n+1}+x_{n-1}$) succesfully. — AtoZ, Feb 21 '19 at 11:52
@Dmitry i am not sure if I could explain to you well, but the code below runs well for the equation without V_n part. Everything else is just definitions of various things which you could ignore. I will be grateful if you could help add the V_n term. You can just write an answer and I will accept. Thanks. — AtoZ, Feb 21 '19 at 11:54

score 0 · Answer 1 · answered Feb 21 '19 at 11:41

function testODE
k=pi/2; p=[-40:40];
TT=1;
z=zeros(1,numel(p)-2);
w=-2*cos(k);
gamma=1;
z=zeros(1,numel(p)-2)
v=-2.5;
v1=-2.625;
v2=-2.375;
psi3=TT*exp(3*i*k);
psi2=TT*exp(2*i*k);
psi1=-psi3+(v2-w+gamma*abs(psi2).^2).*psi2;
psi0=-psi2+(v1-w+gamma*abs(psi1).^2).*psi1;
R0=(psi0.*exp(-i*k)-psi1)./(exp(-i*k)-exp(i*k));
R1=(psi0.*exp(i*k)-psi1)./(exp(i*k)-exp(-i*k));
y = (p <= 1) .* (R0 * exp(1i*k*p) + R1 * exp(-1i*k*p)) + (p >= 2) .* (TT * exp(1i*k*p));
z(1+(numel(p)-2)/2:2+(numel(p)-2)/2)=[v*y(1+ceil(end/2)),v*y(2+ceil(end/2))]
f = @(t,y) [(R0 * exp(1i*k*p(1)) + R1 * exp(-1i*k*p(1))).*exp(-1i*w*t); diff(y,2)+2*y(2:end-1); (TT * exp(1i*k*p(end))).*exp(-1i*w*t)];
y0 = y;
t0 = 0;
T = 5;
N = 50;
[Y,t] = RK42d(f, t0, T, y0, N);
figure;
hold on;
plot(p,Y,'r');
plot(p,Y,'b');
end

function [Y, t] = RK42d(f, t0, T, y0, N)
%the input of f and y0 must both be vectors
%composed of two fuctions and their respective
%intial conditions
%the vector Y will return a vector as well
h = (T - t0)/(N - 1); %Calculate and store the step size
Y = zeros(N,81); %Initialize the X and Y vector
t = linspace(t0,T,N); % A vector to store the time values
Y(1,:) = y0'; % Start Y vector at the intial values.
for i = 1:(N-1)
    y = Y(i,:)';
    k1 = f(t(i),y);
    k2 = f(t(i) +0.5*h, y +0.5*h*k1);
    k3 = f(t(i) +0.5*h, y +0.5*h*k2);
    k4 = f(t(i) +h    , y +h*k3);
    Y(i+1,:) = y + (h/6)*(k1+ 2.*k2 + 2*k3 + k4);
    %Update approximation
end
end

In principle: to add $V_n$ (which is coded as 'z') with $x_{n+1}+x_{n-1}$ (which is coded as diff(y,2)+2*y(2:end-1)' )in "f = @(t,y)......."? — AtoZ, Feb 21 '19 at 12:22

Solving system of ODEs in Matlab

1 Answers1