The differential equation problem with variable solution by using ode45

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

[Walter Roberson]

If you have the symbolic toolbox then I recommend that you look at the flow of calls in the first example under the odeFunction documentation.

[abhishek singh]

how I will get x[t] in that case .

[darova]

Solve this equation first

Use results to solve the system of ODE

[abhishek singh]

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

Answer 1

[darova]

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

[abhishek singh]

This is the matlab file i have constructed what else i have to change in this.noscllator have second order diffrential equation and coupled equation is in xotss. what else i have to do.

[darova]

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

[abhishek singh]

Thank you for helping me.

Answer 2

[abhishek singh]

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

[darova]

What the error says?

Answer 3

[abhishek singh]

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)';

[darova]

It means

सदस्यता सूचकांकों को वास्तविक धनात्मक पूर्णांक या तार्किक होना चाहिए।

In your language. Any ideas what the problem it might be?

Answer 4

[abhishek singh]

I am new here so i don't know what the problem is.

Answer 5

[Walter Roberson]

[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

[darova]

Thanks bro! It works fine!

[abhishek singh]

Thanks

[Walter Roberson]

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) .

Answer 6

[abhishek singh]

No , I am not using derivative ,it complex conjugate only.

[Rik]

Please do not post your comments as answer. Their order can change, which makes it confusing.

Answer 7

[abhishek singh]

I am not getting plot

plot(t2,rho11).

[Walter Roberson]

I wonder: do you happen to be using MacOS Catalina ?