solving differential equation for diffrent intial condition

1 vue (au cours des 30 derniers jours)
samee baig
samee baig le 14 Juil 2016
Commenté : samee baig le 14 Juil 2016
I am trying to solve differential equation for different initial conditions and plotting the result for each initial condition.i need help on writing a code using a loop to do the same and plotting the result against each initial condition. my current code is as below function xdot=fifteenstate2(t,x)
T=1;
v(1) =-x(1)-x(2)*x(3)-x(3)*x(4)-x(5)*x(6)-x(7)*x(8)-x(9)*x(10)-x(11)*x(12)-x(13)*x(14);
v(2) = -x(2);
v(3) =-x(5);
v(4) =-x(5)*x(6)-x(7);
v(5) =-x(9);
v(6) =-x(11);
v(7) =-x(13);
v(8) =x(3);
v(9) =x(6);
v(10) =x(8);
v(11) =x(10);
v(12) =x(12);
v(13) =x(14);
v(14) =x(4);
v=[v(1);v(2);v(3);v(4);v(5);v(6);v(7);v(8);v(9);v(10);v(11);v(12);v(13);v(14)];
%uu=[1;1;1;1;1;1;1];
d(1)=3.54491*sign(v(8))*sqrt(abs(v(8))/T);
d(2)=3.54491*sign(v(9))*sqrt(abs(v(9))/T);
d(3)=3.54491*sign(v(10))*sqrt(abs(v(10))/T);
d(4)=3.54491*sign(v(11))*sqrt(abs(v(11))/T);
d(5)=3.54491*sign(v(12))*sqrt(abs(v(12))/T);
d(6)=3.54491*sign(v(13))*sqrt(abs(v(13))/T);
d(7)=3.54491*sign(v(14))*sqrt(abs(v(14))/T);
d=[d(1);d(2);d(3);d(4);d(5);d(6);d(7)];
y(1)=sin((2*pi*t)/T); y(2)=cos((2*pi*t)/T);
y=[y(1);y(2)];
uu(1)=v(1)+d(1)*y(2)+d(2)*y(2)-d(3)*y(2)-d(4)*y(1)-d(5)*y(1)+d(6)*y(2);%+d(7)*y(1)+d(1)*d(7)*y(2)*y(2);
uu(2)=v(2)+d(1)*y(1);%-d(1)*d(7)*y(1)*y(2);
uu(3)=v(3)+d(2)*y(1);
uu(4)=v(4)+d(3)*y(1);
uu(5)=v(5)+d(4)*y(2);
uu(6)=v(6)+d(5)*y(2);
uu(7)=v(7)+d(6)*y(1);
uu=[uu(1);uu(2);uu(3);uu(4);uu(5);uu(6);uu(7)];
xdot(1)=uu(1); xdot(2)=uu(2); xdot(3)=x(2)*uu(1); xdot(4)=x(3)*uu(1); xdot(5)=uu(3); xdot(6)=x(5)*uu(1); xdot(7)=uu(4); xdot(8)=x(7)*uu(1); xdot(9)=uu(5); xdot(10)=x(9)*uu(1); xdot(11)=uu(6); xdot(12)=x(11)*uu(1); xdot(13)=uu(7); xdot(14)=x(13)*uu(1); xdot=xdot';
clear all close all clc
Tspan =linspace(0,100,500);
IC=[1;-3;2;-2;3;-4;1;-2;5;3;-2;-2;2;-3]; %options = odeset('RelTol',1e-2,'AbsTol',[1e-2 1e-2 1e-2 1e-2 1e-2 1e-2]); [T,Y] = ode45(@fifteenstate2, Tspan,IC)
x1 = Y(:,1); x2 = Y(:,2); x3 = Y(:,3); x4 = Y(:,4); x5 = Y(:,5); x6 = Y(:,6); x7 = Y(:,7); x8 = Y(:,8); x9 = Y(:,9); x10 = Y(:,10); x11 = Y(:,11); x12 = Y(:,12); x13 = Y(:,13); x14 = Y(:,14);
figure(1) plot(T,x1,T,x3,T,x4,T,x6,T,x8,T,x10,T,x12,T,x14,'linewidth',2) legend('x1','x3','x4','x6','x8','x10','x12','x14',8); xlabel('t(s)') ylabel('x_1,x_3,x_4,x_6,x_8,x_10,x_12,x_14') grid on
figure(2) plot(T,x2,T,x5,T,x7,T,x9,T,x11,T,x13,'linewidth',2) legend('x2','x5','x7','x9','x11','x13',6); xlabel('t(s)') ylabel('x_2,x_5,x_7,x_9,x_11,x_13') grid on
  2 commentaires
KSSV
KSSV le 14 Juil 2016
ok.....go ahead, do it.......is some thing objecting you?
samee baig
samee baig le 14 Juil 2016
yes,cannot seem to figure out how to write the code

Connectez-vous pour commenter.

Réponse acceptée

Mischa Kim
Mischa Kim le 14 Juil 2016
Samee, check out this answer. It does exactly what you are looking for.
  1 commentaire
samee baig
samee baig le 14 Juil 2016
its not working,is it because the initial condition i wrote are in vector form ?

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Mathematics dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by