help solve odefcn with time dependence

19 Nov 2022
3 Answers
Answer Accepted
Updated 19 Nov 2022
11 Views (30 days)

0 votes

Hi,
Can you please help me solve my problem with this 16 coupled diffrential equation ? With Jkl-all time dependent series of size (1,991).
Thank you very much!
function dydt = odecdn(t,x,A,U,J,N)
q1= real(angle(conj(A(:,3)).*A(:,1))./(abs(A(:,1)).*sqrt(N-abs(A(:,1)).^2-abs(A(:,2)).^2)));
q2= real(angle(conj(A(:,3)).*A(:,2))./(abs(A(:,1)).*sqrt(N-abs(A(:,2)).^2-abs(A(:,2)).^2)));
p1 = (abs(A(:,1))).^2;
p2 = (abs(A(:,2))).^2;
J11= transpose(real(sqrt(p2./p1).*sin(q2-q1).*(-J/2)));
J12= transpose(real((-J/2)*(-sqrt(p2./p1).*sin(q2-q1)+ sqrt((p2)./(N-p1-p2)).*sin(q2))));
J13= transpose(real((-J/2).*(cos(q2-q1).*sqrt(p2).*(-0.5).*p1.^(-3.2) - cos(q2).*sqrt(p2).*(-0.5).*(N-p1-p2).^(-3/2).*(-1))+2.*U));
J14= transpose(real(0.5.*(p2).^(-0.5).*(p1).^(-0.5).*cos(q2-q1)-(0.5*(p2/(N-p1-p2)).^(-0.5).*cos(q2).*(((N-p1))./((N-p1-p2).^2))) + U));
J21= transpose(real((-J/2).*sin(q2-q1)));
J22= transpose(real((-J/2).*(-sin(q2-q1)+sin(q2).*((N-p1-2.*p2)./(sqrt(p2.*(N-p1-p2)))))));
J23= transpose(real((-J/2).*(-cos(q2).*((-sqrt(p2.*(N-p1-p2))-0.5.*(p2.*(N-p1-p2)).^(-3./2).*(N-p1-2.*p2))./(p2.*(N-p1-p2))))+U));
J24= transpose(real((-J/2).*(-(-2.*sqrt(p2.*(N-p1-p2))-(0.5.*cos(q2).*(N-p1-2.*p2).*(p2.*(N-p1-p2)).^(-3./2))))+2.*U));
J31= transpose(real(-J*sqrt(p1.*p2).*cos(q1-q2)));
J32= transpose(real(J*sqrt(p1.*p2).*cos(q1-q2)));
J33= transpose(real(-J*0.5.*sqrt(p2).*(p1).^(-0.5).*sin(q1-q2)));
J34= transpose(real(-0.5.*J*sqrt(p1).*(p2).^(-0.5).*sin(q1-q2)));
J41= transpose(real(J*sqrt(p1.*p2).*cos(q1-q2)));
J42= transpose(real(J*(-sqrt(p1.*p2).*cos(q1-q2)+sqrt(p2.*(N-p1-p2)).*cos(q2))));
J43= transpose(real(J*(sqrt(p2).*0.5.*p1.^(-0.5).*sin(q1-q2)+sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5))));
J44= transpose(real(J*(sqrt(p1).*0.5.*(p2).^(-0.5).*sin(q1-q2)-sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5).*(-p1-2.*p2))));
dydt=zeros(9991,9991)
y=zeros(16,9991)
dydt(1)= (J11.*(y(1))+ (J12).*(y(5))+(J13).*(y(9))+(J14).*(y(13)));
dydt(2)= J11.*(y(2))+ (J12).*(y(6))+(J13).*(y(10))+(J14).*(y(14));
dydt(3)= (J11.*(y(3))+ J12.*(y(7))+J13.*(y(11))+J14.*(y(15)));
dydt(4)= (J11.*(y(4))+ J12.*(y(8))+J13.*(y(12))+J14.*(y(16)));
dydt(5)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
dydt(6)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
dydt(7)= (J21.*(y(3))+ J22.*(y(7))+J23.*(y(11))+J24.*(y(15)));
dydt(8)= (J21.*(y(4))+ J22.*(y(8))+J23.*(y(12))+J24.*(y(16)));
dydt(9)= (J31.*(y(1))+ J32.*(y(5))+J33.*(y(9))+J34.*(y(13)));
dydt(10)= (J31.*(y(2))+ J32.*(y(6))+J33.*(y(10))+J34.*(y(14)));
dydt(11)= (J31.*(y(3))+ J32.*(y(7))+J33.*(y(11))+J34.*(y(15)));
dydt(12)= (J31.*(y(4))+ J32.*(y(8))+J33.*(y(12))+J34.*(y(16)));
dydt(13)= (J41.*(y(1))+ J42.*(y(5))+J43.*(y(9))+J44.*(y(13)));
dydt(14)= (J41.*(y(2))+ J42.*(y(6))+J43.*(y(10))+J44.*(y(14)));
dydt(15)= (J41.*(y(3))+ J42.*(y(7))+J43.*(y(11))+J44.*(y(15)));
dydt(16)= (J41.*(y(4))+ J42.*(y(8))+J43.*(y(12))+J44.*(y(16)));
end

1 Comment
Show -1 older comments Hide -1 older comments

Torsten
Torsten on 19 Nov 2022
The number of differential equations and the number of unknown functions must be the same.
Since you initialize
dydt=zeros(9991,9991)
y=zeros(16,9991)
you seem to have 16*9991 unknown functions and 9991*9991 differential equations.
This is not solvable.

Sign in to comment.

Sign in to answer this question.

 Accepted Answer

Torsten
Torsten on 19 Nov 2022

0 votes

tspan = 1:0.1:10^3;
y0 = [sqrt(0.9) sqrt(1.1) sqrt(1)];
U=1;
J=1;
N=3; % gewenste aantal deeltjes per site
mu=U*N-J; % chemische potentiaal
%lambda=J/UN
%k=J
%abs(y(length(t),1))^2+abs(y(length(t),2))^2
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
% grondtoestand vinden
[t,y] = ode15s(@(t,y) odefcnim(y,U,J,mu,N), tspan, y0, options);
%plot(t,abs(y(:,1:3)).^2)
% begintoestand
y0=y(end,:).*exp(1i*[1 0.5 -1].*0.1.*[sqrt(0.9) sqrt(1.1) sqrt(1)]);
% dynamica
[t,y] = ode15s(@(t,y) odefcn(y,U,J,N), tspan, y0, options);
%solve for A such that, dA/dt= JA , with J the jacobian matrix.
A=y;
%y=A=9991x3-matrix
tspan = [1:0.1:10^3];
ic = [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ];
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
[t,s] = ode15s(@(t,y) odefcn2(t,y,tspan,A,U,J,N), tspan, ic, opts);
%odesolver calculates A in function of the time.
function dydt = odefcn(y,U,J,N)
dydt = zeros(3,1);
dydt(1) = -1i * (-J*0.5*y(2)+ U*N*abs(y(1))^2*y(1));
dydt(2) = -1i * (-J*0.5*y(1)-J*0.5*y(3)+ U*N*abs(y(2))^2*y(2));
dydt(3) = -1i * (-J*0.5*y(2)+U*N*(N-norm(y(1))^2-norm(y(2))^2)*y(3));
end
function dydt = odefcnim(y,U,J,mu,N)
dydt = zeros(3,1);
dydt(1) = -(-J*0.5*y(2)+ U*N*(abs(y(1))^2)*y(1) - mu*y(1) );
dydt(2) = -(-J*0.5*y(1)-J*0.5*y(3)+ U*N*(abs(y(2))^2)*y(2) - mu*y(2) );
dydt(3) = -(-J*0.5*y(2)+U*N*((N-norm(y(1))^2-norm(y(2))^2))*y(3) - mu*y(3) );
end
function dydt = odefcn2(t,y,tspan,Amat,U,J,N)
A = interp1(tspan,Amat,t);
q1= real(angle(conj(A(:,3)).*A(:,1))./(abs(A(:,1)).*sqrt(N-abs(A(:,1)).^2-abs(A(:,2)).^2)));
q2= real(angle(conj(A(:,3)).*A(:,2))./(abs(A(:,1)).*sqrt(N-abs(A(:,2)).^2-abs(A(:,2)).^2)));
p1 = (abs(A(:,1))).^2;
p2 = (abs(A(:,2))).^2;
J11= transpose(real(sqrt(p2./p1).*sin(q2-q1).*(-J/2)));
J12= transpose(real((-J/2)*(-sqrt(p2./p1).*sin(q2-q1)+ sqrt((p2)./(N-p1-p2)).*sin(q2))));
J13= transpose(real((-J/2).*(cos(q2-q1).*sqrt(p2).*(-0.5).*p1.^(-3.2) - cos(q2).*sqrt(p2).*(-0.5).*(N-p1-p2).^(-3/2).*(-1))+2.*U));
J14= transpose(real(0.5.*(p2).^(-0.5).*(p1).^(-0.5).*cos(q2-q1)-(0.5*(p2/(N-p1-p2)).^(-0.5).*cos(q2).*(((N-p1))./((N-p1-p2).^2))) + U));
J21= transpose(real((-J/2).*sin(q2-q1)));
J22= transpose(real((-J/2).*(-sin(q2-q1)+sin(q2).*((N-p1-2.*p2)./(sqrt(p2.*(N-p1-p2)))))));
J23= transpose(real((-J/2).*(-cos(q2).*((-sqrt(p2.*(N-p1-p2))-0.5.*(p2.*(N-p1-p2)).^(-3./2).*(N-p1-2.*p2))./(p2.*(N-p1-p2))))+U));
J24= transpose(real((-J/2).*(-(-2.*sqrt(p2.*(N-p1-p2))-(0.5.*cos(q2).*(N-p1-2.*p2).*(p2.*(N-p1-p2)).^(-3./2))))+2.*U));
J31= transpose(real(-J*sqrt(p1.*p2).*cos(q1-q2)));
J32= transpose(real(J*sqrt(p1.*p2).*cos(q1-q2)));
J33= transpose(real(-J*0.5.*sqrt(p2).*(p1).^(-0.5).*sin(q1-q2)));
J34= transpose(real(-0.5.*J*sqrt(p1).*(p2).^(-0.5).*sin(q1-q2)));
J41= transpose(real(J.*sqrt(p1.*p2).*cos(q1-q2)));
J42= transpose(real(J.*(-sqrt(p1.*p2).*cos(q1-q2)+sqrt(p2.*(N-p1-p2)).*cos(q2))));
J43= transpose(real(J.*(sqrt(p2).*0.5.*p1.^(-0.5).*sin(q1-q2)+sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5))));
J44= transpose(real(J.*(sqrt(p1).*0.5.*(p2).^(-0.5).*sin(q1-q2)-sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5).*(-p1-2.*p2))));
dydt = zeros(16,1);
dydt(1)= (J11.*(y(1))+ (J12).*(y(5))+(J13).*(y(9))+(J14).*(y(13)));
dydt(2)= J11.*(y(2))+ (J12).*(y(6))+(J13).*(y(10))+(J14).*(y(14));
dydt(3)= (J11.*(y(3))+ J12.*(y(7))+J13.*(y(11))+J14.*(y(15)));
dydt(4)= (J11.*(y(4))+ J12.*(y(8))+J13.*(y(12))+J14.*(y(16)));
dydt(5)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
dydt(6)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
dydt(7)= (J21.*(y(3))+ J22.*(y(7))+J23.*(y(11))+J24.*(y(15)));
dydt(8)= (J21.*(y(4))+ J22.*(y(8))+J23.*(y(12))+J24.*(y(16)));
dydt(9)= (J31.*(y(1))+ J32.*(y(5))+J33.*(y(9))+J34.*(y(13)));
dydt(10)= (J31.*(y(2))+ J32.*(y(6))+J33.*(y(10))+J34.*(y(14)));
dydt(11)= (J31.*(y(3))+ J32.*(y(7))+J33.*(y(11))+J34.*(y(15)));
dydt(12)= (J31.*(y(4))+ J32.*(y(8))+J33.*(y(12))+J34.*(y(16)));
dydt(13)= (J41.*(y(1))+ J42.*(y(5))+J43.*(y(9))+J44.*(y(13)));
dydt(14)= (J41.*(y(2))+ J42.*(y(6))+J43.*(y(10))+J44.*(y(14)));
dydt(15)= (J41.*(y(3))+ J42.*(y(7))+J43.*(y(11))+J44.*(y(15)));
dydt(16)= (J41.*(y(4))+ J42.*(y(8))+J43.*(y(12))+J44.*(y(16)));
end

More Answers (2)

Hélène De Witte
Hélène De Witte on 19 Nov 2022
Edited: Hélène De Witte on 19 Nov 2022

0 votes

Hi, I have tried dydt=zeros(9991,9991) y=zeros(16,9991) , tol, but even then it blocked. But the Jkl's are time dependent is there an option to get a counter,I have to solve: dM/dt = J*M and J and M are 4x4 matrices, and J is time fependent,I only know matrix J in function of the time-step. So I wrote the 16 deffrential equations but when I run it it blocks. Can you please help me a second time, please?
I'm really thankful! Thank you very much!
%solve for A such that, dA/dt= JA , with J the jacobian matrix.
A=y;
%y=A=9991x3-matrix
tspan = [1:0.1:10^3];
ic = [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ];
opts = odeset('RelTol',1e-20,'AbsTol',1e-20);
[t,s] = ode45(@(t,y) odefcn(t,y,A,U,J,N), tspan, ic, opts);
%odesolver calculates A in function of the time.
function dydt = odefcn(t,x,A,U,J,N)
q1= real(angle(conj(A(:,3)).*A(:,1))./(abs(A(:,1)).*sqrt(N-abs(A(:,1)).^2-abs(A(:,2)).^2)));
q2= real(angle(conj(A(:,3)).*A(:,2))./(abs(A(:,1)).*sqrt(N-abs(A(:,2)).^2-abs(A(:,2)).^2)));
p1 = (abs(A(:,1))).^2;
p2 = (abs(A(:,2))).^2;
J11= transpose(real(sqrt(p2./p1).*sin(q2-q1).*(-J/2)));
J12= transpose(real((-J/2)*(-sqrt(p2./p1).*sin(q2-q1)+ sqrt((p2)./(N-p1-p2)).*sin(q2))));
J13= transpose(real((-J/2).*(cos(q2-q1).*sqrt(p2).*(-0.5).*p1.^(-3.2) - cos(q2).*sqrt(p2).*(-0.5).*(N-p1-p2).^(-3/2).*(-1))+2.*U));
J14= transpose(real(0.5.*(p2).^(-0.5).*(p1).^(-0.5).*cos(q2-q1)-(0.5*(p2/(N-p1-p2)).^(-0.5).*cos(q2).*(((N-p1))./((N-p1-p2).^2))) + U));
J21= transpose(real((-J/2).*sin(q2-q1)));
J22= transpose(real((-J/2).*(-sin(q2-q1)+sin(q2).*((N-p1-2.*p2)./(sqrt(p2.*(N-p1-p2)))))));
J23= transpose(real((-J/2).*(-cos(q2).*((-sqrt(p2.*(N-p1-p2))-0.5.*(p2.*(N-p1-p2)).^(-3./2).*(N-p1-2.*p2))./(p2.*(N-p1-p2))))+U));
J24= transpose(real((-J/2).*(-(-2.*sqrt(p2.*(N-p1-p2))-(0.5.*cos(q2).*(N-p1-2.*p2).*(p2.*(N-p1-p2)).^(-3./2))))+2.*U));
J31= transpose(real(-J*sqrt(p1.*p2).*cos(q1-q2)));
J32= transpose(real(J*sqrt(p1.*p2).*cos(q1-q2)));
J33= transpose(real(-J*0.5.*sqrt(p2).*(p1).^(-0.5).*sin(q1-q2)));
J34= transpose(real(-0.5.*J*sqrt(p1).*(p2).^(-0.5).*sin(q1-q2)));
J41= transpose(real(J.*sqrt(p1.*p2).*cos(q1-q2)));
J42= transpose(real(J.*(-sqrt(p1.*p2).*cos(q1-q2)+sqrt(p2.*(N-p1-p2)).*cos(q2))));
J43= transpose(real(J.*(sqrt(p2).*0.5.*p1.^(-0.5).*sin(q1-q2)+sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5))));
J44= transpose(real(J.*(sqrt(p1).*0.5.*(p2).^(-0.5).*sin(q1-q2)-sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5).*(-p1-2.*p2))));
dydt=zeros(16,9991)
y=zeros(16,9991)
dydt(1)= (J11.*(y(1))+ (J12).*(y(5))+(J13).*(y(9))+(J14).*(y(13)));
dydt(2)= J11.*(y(2))+ (J12).*(y(6))+(J13).*(y(10))+(J14).*(y(14));
dydt(3)= (J11.*(y(3))+ J12.*(y(7))+J13.*(y(11))+J14.*(y(15)));
dydt(4)= (J11.*(y(4))+ J12.*(y(8))+J13.*(y(12))+J14.*(y(16)));
dydt(5)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
dydt(6)= (J21.*(y(1))+ J22.*(y(5))+J23.*(y(9))+J24.*(y(13)));
dydt(7)= (J21.*(y(3))+ J22.*(y(7))+J23.*(y(11))+J24.*(y(15)));
dydt(8)= (J21.*(y(4))+ J22.*(y(8))+J23.*(y(12))+J24.*(y(16)));
dydt(9)= (J31.*(y(1))+ J32.*(y(5))+J33.*(y(9))+J34.*(y(13)));
dydt(10)= (J31.*(y(2))+ J32.*(y(6))+J33.*(y(10))+J34.*(y(14)));
dydt(11)= (J31.*(y(3))+ J32.*(y(7))+J33.*(y(11))+J34.*(y(15)));
dydt(12)= (J31.*(y(4))+ J32.*(y(8))+J33.*(y(12))+J34.*(y(16)));
dydt(13)= (J41.*(y(1))+ J42.*(y(5))+J43.*(y(9))+J44.*(y(13)));
dydt(14)= (J41.*(y(2))+ J42.*(y(6))+J43.*(y(10))+J44.*(y(14)));
dydt(15)= (J41.*(y(3))+ J42.*(y(7))+J43.*(y(11))+J44.*(y(15)));
dydt(16)= (J41.*(y(4))+ J42.*(y(8))+J43.*(y(12))+J44.*(y(16)));
end

2 Comments
Show None Hide None

Torsten
Torsten on 19 Nov 2022
Edited: Torsten on 19 Nov 2022
What are the time instants that correspond to the 9991 entries of your matrix A ?
Hélène De Witte
Hélène De Witte on 19 Nov 2022
Thank you very much! Well, it's also a solution (I try to solve a biger problem and needed the y first, it were also diffrential equations) of a ode function defined as follow:
tspan = 1:0.1:10^3;
y0 = [sqrt(0.9) sqrt(1.1) sqrt(1)];
U=1;
J=1;
N=3; % gewenste aantal deeltjes per site
mu=U*N-J; % chemische potentiaal
%lambda=J/UN
%k=J
%abs(y(length(t),1))^2+abs(y(length(t),2))^2
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
% grondtoestand vinden
[t,y] = ode15s(@(t,y) odefcnim(y,U,J,mu,N), tspan, y0, options);
%plot(t,abs(y(:,1:3)).^2)
% begintoestand
y0=y(end,:).*exp(1i*[1 0.5 -1].*0.1.*[sqrt(0.9) sqrt(1.1) sqrt(1)]);
% dynamica
[t,y] = ode15s(@(t,y) odefcn(y,U,J,N), tspan, y0, options);
function dydt = odefcn(y,U,J,N)
dydt = zeros(3,1);
dydt(1) = -1i * (-J*0.5*y(2)+ U*N*abs(y(1))^2*y(1));
dydt(2) = -1i * (-J*0.5*y(1)-J*0.5*y(3)+ U*N*abs(y(2))^2*y(2));
dydt(3) = -1i * (-J*0.5*y(2)+U*N*(N-norm(y(1))^2-norm(y(2))^2)*y(3));
end
function dydt = odefcnim(y,U,J,mu,N)
dydt = zeros(3,1);
dydt(1) = -(-J*0.5*y(2)+ U*N*(abs(y(1))^2)*y(1) - mu*y(1) );
dydt(2) = -(-J*0.5*y(1)-J*0.5*y(3)+ U*N*(abs(y(2))^2)*y(2) - mu*y(2) );
dydt(3) = -(-J*0.5*y(2)+U*N*((N-norm(y(1))^2-norm(y(2))^2))*y(3) - mu*y(3) );
end

Sign in to comment.

Hélène De Witte
Hélène De Witte on 19 Nov 2022

0 votes

Thank you very much! Maybe this helps too, In fact, I can solve it, like this, but I want to understand how I can fix it without manupulating the problem, because it ask a lot of time...
s=[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]
for stap=1:2
teller=1
tspan = [teller:0.1:teller+0.1];
ic = [s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end) s(end)];
opts = odeset('RelTol',1e-20,'AbsTol',1e-20);
[t,s] = ode45(@(t,y) kale(t,y,A,U,J,N,stap), tspan, ic, opts);
%odesolver calculates A in function of the time.
teller=teller+0.1
end
function dydt = kale(t,y,A,U,J,N, stap)
q1= real(angle(conj(A(:,3)).*A(:,1))./(abs(A(:,1)).*sqrt(N-abs(A(:,1)).^2-abs(A(:,2)).^2)));
q2= real(angle(conj(A(:,3)).*A(:,2))./(abs(A(:,1)).*sqrt(N-abs(A(:,2)).^2-abs(A(:,2)).^2)));
p1 = (abs(A(:,1))).^2;
p2 = (abs(A(:,2))).^2;
J11= (real(sqrt(p2./p1).*sin(q2-q1).*(-J/2)));
J12= (real((-J/2)*(-sqrt(p2./p1).*sin(q2-q1)+ sqrt((p2)./(N-p1-p2)).*sin(q2))));
J13= (real((-J/2).*(cos(q2-q1).*sqrt(p2).*(-0.5).*p1.^(-3.2) - cos(q2).*sqrt(p2).*(-0.5).*(N-p1-p2).^(-3/2).*(-1))+2.*U));
J14= (real(0.5.*(p2).^(-0.5).*(p1).^(-0.5).*cos(q2-q1)-(0.5*(p2/(N-p1-p2)).^(-0.5).*cos(q2).*(((N-p1))./((N-p1-p2).^2))) + U));
J21= (real((-J/2).*sin(q2-q1)));
J22= (real((-J/2).*(-sin(q2-q1)+sin(q2).*((N-p1-2.*p2)./(sqrt(p2.*(N-p1-p2)))))));
J23= (real((-J/2).*(-cos(q2).*((-sqrt(p2.*(N-p1-p2))-0.5.*(p2.*(N-p1-p2)).^(-3./2).*(N-p1-2.*p2))./(p2.*(N-p1-p2))))+U));
J24= (real((-J/2).*(-(-2.*sqrt(p2.*(N-p1-p2))-(0.5.*cos(q2).*(N-p1-2.*p2).*(p2.*(N-p1-p2)).^(-3./2))))+2.*U));
J31= (real(-J*sqrt(p1.*p2).*cos(q1-q2)));
J32= (real(J.*sqrt(p1.*p2).*cos(q1-q2)));
J33= (real(-J*0.5.*sqrt(p2).*(p1).^(-0.5).*sin(q1-q2)));
J34= (real(-0.5.*J.*sqrt(p1).*(p2).^(-0.5).*sin(q1-q2)));
J41= (real(J.*sqrt(p1.*p2).*cos(q1-q2)));
J42= (real(J.*(-sqrt(p1.*p2).*cos(q1-q2)+sqrt(p2.*(N-p1-p2)).*cos(q2))));
J43= (real(J.*(sqrt(p2).*0.5.*p1.^(-0.5).*sin(q1-q2)+sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5))));
J44= (real(J.*(sqrt(p1).*0.5.*(p2).^(-0.5).*sin(q1-q2)-sin(q2).*0.5.*(p2.*(N-p1-p2)).^(-0.5).*(-p1-2.*p2))));
dydt= zeros(16,1);
dydt(1)= (J11(stap,1)*(y(1))+ (J12(stap,1))*(y(5))+(J13(stap,1))*(y(9))+(J14(stap,1))*(y(13)));
dydt(2)= J11(stap,1)*(y(2))+ (J12(stap,1))*(y(6))+(J13(stap,1))*(y(10))+(J14(stap,1))*(y(14));
dydt(3)= (J11(stap,1)*(y(3))+ J12(stap,1)*(y(7))+J13(stap,1)*(y(11))+J14(stap,1)*(y(15)));
dydt(4)= (J11(stap,1)*(y(4))+ J12(stap,1)*(y(8))+J13(stap,1)*(y(12))+J14(stap,1)*(y(16)));
dydt(5)= (J21(stap,1)*(y(1))+ J22(stap,1)*(y(5))+J23(stap,1)*(y(9))+J24(stap,1)*(y(13)));
dydt(6)= (J21(stap,1)*(y(1))+ J22(stap,1)*(y(5))+J23(stap,1)*(y(9))+J24(stap,1)*(y(13)));
dydt(7)= (J21(stap,1)*(y(3))+ J22(stap,1)*(y(7))+J23(stap,1)*(y(11))+J24(stap,1)*(y(15)));
dydt(8)= (J21(stap,1)*(y(4))+ J22(stap,1)*(y(8))+J23(stap,1)*(y(12))+J24(stap,1)*(y(16)));
dydt(9)= (J31(stap,1)*(y(1))+ J32(stap,1)*(y(5))+J33(stap,1)*(y(9))+J34(stap,1)*(y(13)));
dydt(10)= (J31(stap,1)*(y(2))+ J32(stap,1)*(y(6))+J33(stap,1)*(y(10))+J34(stap,1)*(y(14)));
dydt(11)= (J31(stap,1)*(y(3))+ J32(stap,1)*(y(7))+J33(stap,1)*(y(11))+J34(stap,1)*(y(15)));
dydt(12)= (J31(stap,1)*(y(4))+ J32(stap,1)*(y(8))+J33(stap,1)*(y(12))+J34(stap,1)*(y(16)));
dydt(13)= (J41(stap,1)*(y(1))+ J42(stap,1)*(y(5))+J43(stap,1)*(y(9))+J44(stap,1)*(y(13)));
dydt(14)= (J41(stap,1)*(y(2))+ J42(stap,1)*(y(6))+J43(stap,1)*(y(10))+J44(stap,1)*(y(14)));
dydt(15)= (J41(stap,1)*(y(3))+ J42(stap,1)*(y(7))+J43(stap,1)*(y(11))+J44(stap,1)*(y(15)));
dydt(16)= (J41(stap,1)*(y(4))+ J42(stap,1)*(y(8))+J43(stap,1)*(y(12))+J44(stap,1)*(y(16)));
end

Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!