solving differential equation for diffrent intial condition

조회 수: 1 (최근 30일)
samee baig
samee baig 2016년 7월 14일
댓글: samee baig 2016년 7월 14일
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
KSSV
KSSV 2016년 7월 14일
ok.....go ahead, do it.......is some thing objecting you?
samee baig
samee baig 2016년 7월 14일
yes,cannot seem to figure out how to write the code

댓글을 달려면 로그인하십시오.

채택된 답변

Mischa Kim
Mischa Kim 2016년 7월 14일
Samee, check out this answer. It does exactly what you are looking for.
  댓글 수: 1
samee baig
samee baig 2016년 7월 14일
its not working,is it because the initial condition i wrote are in vector form ?

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Mathematics에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by