The differential equation problem with variable solution by using ode45

Asked by abhishek singh on 31 Oct 2019

Latest activity Commented on by Walter Roberson on 1 Nov 2019

Accepted Answer by darova

7 views (last 30 days)

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 Comments

Show 1 older comment

abhishek singh on 31 Oct 2019

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

darova on 31 Oct 2019

Solve this equation first

Use results to solve the system of ODE

abhishek singh on 31 Oct 2019

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

Products

MATLAB

7 Answers

Link

Answer by darova on 31 Oct 2019

Accepted Answer

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 Comments

abhishek singh on 31 Oct 2019

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 on 31 Oct 2019

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 on 31 Oct 2019

Thank you for helping me.

Link

Answer by abhishek singh on 31 Oct 2019

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

1 Comment

darova on 31 Oct 2019

What the error says?

Link

Answer by abhishek singh on 31 Oct 2019

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 Comment

darova on 31 Oct 2019

It means

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

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

Link

Answer by abhishek singh on 31 Oct 2019

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

0 Comments

Link

Answer by Walter Roberson on 31 Oct 2019

[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 Comments

darova on 31 Oct 2019

Thanks bro! It works fine!

abhishek singh on 1 Nov 2019

Thanks

Walter Roberson on 1 Nov 2019

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