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 commentaire
Afficher -1 commentaires plus anciens Masquer -1 commentaires plus anciens

Torsten
Torsten le 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.

Connectez-vous pour commenter.

Connectez-vous pour répondre à cette question.

 Réponse acceptée

Torsten
Torsten le 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

Plus de réponses (2)

Hélène De Witte
Hélène De Witte le 19 Nov 2022
Modifié(e) : Hélène De Witte le 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 commentaires
Afficher Aucune Masquer Aucune

Torsten
Torsten le 19 Nov 2022
Modifié(e) : Torsten le 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 le 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

Connectez-vous pour commenter.

Hélène De Witte
Hélène De Witte le 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

Catégories

En savoir plus sur Ordinary Differential Equations dans Centre d'aide et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by