solving differential equation for diffrent intial condition

Question

0 个投票

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