The differential equation problem with variable solution by using ode45

1 次查看(过去 30 天)
I have four coupled diffrential equation shown bellow :-
In which a ,b ,c ,d ,e,f are constant and x[t] we will get from the solution of this second order diffrential equation .how to write code for it in matlab.plz help
  4 个评论
abhishek singh
abhishek singh 2019-10-31

I have solved individually the second order differential equation how to use it output as x(t) in coupled differential equation.

请先登录,再进行评论。

采纳的回答

darova
darova 2019-10-31
Here is an idea
[t1,x1] = ode45(@F1,ts,x0); % solve second order
[t2,s] = ode45(@(t,s)F2(t,s,t1,x1(:,1)),ts,s0); % solve system
function ds = F2(t,s,t1,x1)
x = interp1(t1,x1,t); % extract x
ds(1,1) = a*x-b ...
ds(2,1) = c*x*s(1) ...
end
  3 个评论
darova
darova 2019-10-31
编辑:darova 2019-10-31
Here is how your code should look like
function main
%% your constants
[t1,x1] = ode45(@noscillator,[0 10],[0 1]);
[t2,s1] = ode45(@(t,s) xotss(t,s,t1,x1(:,1)), [0 10], [1 0 0 0]);
plot(t2,s1(:,1))
function xdot=noscillator(t,x)
xdot(1) = x(2);
xdot(2) = -(omega^2)*x(1)-3*(gamma/m)*x(1)^2 - 4*(beta/m)*x(1)^3 + (V/m)*cos(w *t);
xdot=xdot';
end
function dxdt = xotss(t,s,t1,x1)
x = interp1(t1,x1,t);
dxdt(1) = (a*x-b)*s(1) + c*x*s(2);
% ...
dxdt = -1i*dxdt'
end
end

请先登录,再进行评论。

更多回答(6 个)

abhishek singh
abhishek singh 2019-10-31
Sir ,I need to calculate rho11 at diffrent time but i am getting an error plz help.
function F1
%% your constants
[t1,x1] = ode45(@noscillator,[0 100],[0 1]);
[t2,s1] = ode45(@(t,s) xotss(t,s,t1,x1(:,1)), [0 100], [1 0 0 0]);
for ti = 0:1:100
rho11(ti)=s1(ti,1).*s1(ti,1)'-s1(ti,3).*s1(ti,3)';
rho12(ti)=s1(ti,1).*s1(ti,2)'+s1(ti,3).*s1(ti,4)';
rho21(ti)=s1(ti,2).*s1(ti,1)'+s1(ti,4).*s1(ti,3)';
rho22(ti)=s1(ti,2).*s1(ti,2)'-s1(ti,4).*s1(ti,4)';
end
s1(:,1)
plot(t2,rho11+rho22)
function xdot=noscillator(t,x)
m=1; omega=1; beta=0.01;gamma=0.01; V=.5; w=(3.0)/(2*pi);
xdot(1) = x(2);
xdot(2) = -(omega^2)*x(1)-3*(gamma/m)*x(1)^2 - 4*(beta/m)*x(1)^3 + (V/m)*cos(w *t);
xdot=xdot';
end
function dsdt = xotss(t,s,t1,x1)
x = interp1(t1,x1,t);
dsdt = zeros(4,1); % this creates a empty coloumn
%vector that you can fill with four derivatives:
a=1;b=0.05625;c=0.010825;d=0.5;e=0.5;f=0.01875;h=0.00625;% define constants
dsdt(1) = -1i* s(1) *(a+h*x-b) -1i* c*x* s(2);
dsdt(2) = -1i*c*x *s(1) -1i* (d-e-h*x)*s(2) +1i* f*s(3);
dsdt(3) = 1i* f*s(2) -1i*(e+h*x-d) *s(3)-1i* c*x* s(4);
dsdt(4) = -1i* c*x* s(3) -1i* (b-a-h*x)* s(4);
end
end

abhishek singh
abhishek singh 2019-10-31
Subscript indices must either be real positive integers or logicals.
Error in F1 (line 6)
rho11(ti)=s1(ti,1).*s1(ti,1)'-s1(ti,3).*s1(ti,3)';
  1 个评论
darova
darova 2019-10-31
It means
सदस्यता सूचकांकों को वास्तविक धनात्मक पूर्णांक या तार्किक होना चाहिए।
In your language. Any ideas what the problem it might be?

请先登录,再进行评论。


abhishek singh
abhishek singh 2019-10-31
I am new here so i don't know what the problem is.

Walter Roberson
Walter Roberson 2019-10-31
[t1,x1] = ode45(@noscillator,[0:100],[0 1]);
[t2,s1] = ode45(@(t,s) xotss(t,s,t1,x1(:,1)), [0:100], [1 0 0 0]);
for ti = 0:1:100
rho11(ti+1)=s1(ti+1,1).*s1(ti+1,1)'-s1(ti+1,3).*s1(ti+1,3)';
rho12(ti+1)=s1(ti+1,1).*s1(ti+1,2)'+s1(ti+1,3).*s1(ti+1,4)';
rho21(ti+1)=s1(ti+1,2).*s1(ti+1,1)'+s1(ti+1,4).*s1(ti+1,3)';
rho22(ti+1)=s1(ti+1,2).*s1(ti+1,2)'-s1(ti+1,4).*s1(ti+1,4)';
end
  3 个评论
Walter Roberson
Walter Roberson 2019-11-1
Note that s1(ti+1,1)' means the conjugate complex transpose of s1(ti+1,1) . It is, however, a scalar, so transpose does not make any change. The You are also expecting real-valued results, so the conjugate is probably not makeing any changes. I suspect you are doing the equivalent of squaring the value.
I worry that you might have that that s1(ti+1,1)' is the derivative of s1(ti+1,1) .

请先登录,再进行评论。


abhishek singh
abhishek singh 2019-11-1
No , I am not using derivative ,it complex conjugate only.

abhishek singh
abhishek singh 2019-11-1
I am not getting plot
plot(t2,rho11).

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by