plot from ode45

2 views (last 30 days)
shiv gaur
shiv gaur on 3 Feb 2022
Commented: shiv gaur on 3 Feb 2022
function kk1
clear all
close all
options = odeset('RelTol',1e-9,'AbsTol', 1e-16);
a0=0.1;
b0=0.1;
p=0.2;
theta=30:70;
alpha=sqrt(1-0.25*p^2);
a=(alpha-sqrt(1-(p)^2))/sqrt(alpha^2-0.25*sin(theta));
sol =ode45(@q,[0 1200],[0.01 0.01 0.01 0.01 0.01 0], options)
end
% output
x1=sol(1).y(1,:);
y1=sol(1).y(2,:);
plot(x,y)
function dy=q(x)
f=sqrt(1+2*x);
b=(1+2*x)^-1.5;
m=16*alpha^2-0.50*sin(theta);
m=4*(2*alpha^2-8*(alpha*sqrt(1-p^2)));
o=16*sqrt(alpha^2-0.50*sin(theta));
A=y*(1+2.*x)^-1.5.*(2+a);
B=-1i*4*y*(0.633)^2*(alpha.^2-0.25*sin(theta))+2*sqrt(alpha^.2-0.25.*sin(theta))*(alpha-sqrt(1-(p).^2));
C=y*0.633^2*(m+n+o).*(alpha-2*sqrt(1-(p).^2));
D=- y*((2+a)^-1)*f^-2*(1.5+1i*(9/8)*f*b);
E=-0.25*0.633*0.633*y*(4-5*p^2+4*(1-exp(0.18/4)));
F=(1/8)*p^2*b0*a0*((1.70*sin(theta)+3.5*sqrt(alpha^2-0.25*sin(theta))+2*(alpha-sqrt(1-(p)^2))));
dy=A+B+C+D+E+F;
end
end
pl you are req to plot x vs y by varing theta from 30 to 70

Answers (1)

Yongjian Feng
Yongjian Feng on 3 Feb 2022
Too many errors in your code. I tried to clean some, but it still has trouble to figure out what y is:
function kk1
clear all
close all
options = odeset('RelTol',1e-9,'AbsTol', 1e-16);
a0=0.1;
b0=0.1;
p=0.2;
theta=30:70;
alpha=sqrt(1-0.25*p^2);
a=(alpha-sqrt(1-(p)^2))./sqrt(alpha^2-0.25*sin(theta));
sol =ode45(@q,[0 1200],[0.01 0.01 0.01 0.01 0.01 0], options)
x1=sol(1).y(1,:);
y1=sol(1).y(2,:);
plot(x,y)
end
% % output
% x1=sol(1).y(1,:);
% y1=sol(1).y(2,:);
% plot(x,y)
%
function dy=q(x, alpha)
f=sqrt(1+2*x);
b=(1+2*x)^-1.5;
theta=30:70;
p=0.2;
m=16*alpha.^2-0.50*sin(theta);
m=4*(2*alpha.^2-8*(alpha*sqrt(1-p^2)));
o=16*sqrt(alpha.^2-0.50*sin(theta));
A=y*(1+2.*x)^-1.5.*(2+a); %What is y here?
B=-1i*4*y*(0.633)^2*(alpha.^2-0.25*sin(theta))+2*sqrt(alpha^.2-0.25.*sin(theta))*(alpha-sqrt(1-(p).^2));
C=y*0.633^2*(m+n+o).*(alpha-2*sqrt(1-(p).^2));
D=- y*((2+a)^-1)*f^-2*(1.5+1i*(9/8)*f*b);
E=-0.25*0.633*0.633*y*(4-5*p^2+4*(1-exp(0.18/4)));
F=(1/8)*p^2*b0*a0*((1.70*sin(theta)+3.5*sqrt(alpha^2-0.25*sin(theta))+2*(alpha-sqrt(1-(p)^2))));
dy=A+B+C+D+E+F;
end
% end
  2 Comments
shiv gaur
shiv gaur on 3 Feb 2022
edit code
function kk1
clear all
close all
options = odeset('RelTol',1e-9,'AbsTol', 1e-16);
a0=0.1;
b0=0.1;
p=0.2;
theta=30:70;
alpha=sqrt(1-0.25*p^2);
a=(alpha-sqrt(1-(p)^2))./sqrt(alpha^2-0.25*sin(theta));
sol =ode45(@q,[0 1200],0.01] options)
x1=sol(1).y(1,:);
y1=sol(1).y(2,:);
plot(x,y)
end
% % output
% x1=sol(1).y(1,:);
% y1=sol(1).y(2,:);
% plot(x,y)
%
function dy=q(x, alpha)
f=sqrt(1+2*x);
b=(1+2*x)^-1.5;
theta=30:70;
p=0.2;
m=16*alpha.^2-0.50*sin(theta);
m=4*(2*alpha.^2-8*(alpha*sqrt(1-p^2)));
o=16*sqrt(alpha.^2-0.50*sin(theta));
A=y*(1+2.*x)^-1.5.*(2+a); %What is y here?
B=-1i*4*y*(0.633)^2*(alpha.^2-0.25*sin(theta))+2*sqrt(alpha^.2-0.25.*sin(theta))*(alpha-sqrt(1-(p).^2));
C=y*0.633^2*(m+n+o).*(alpha-2*sqrt(1-(p).^2));
D=- y*((2+a)^-1)*f^-2*(1.5+1i*(9/8)*f*b);
E=-0.25*0.633*0.633*y*(4-5*p^2+4*(1-exp(0.18/4)));
F=(1/8)*p^2*b0*a0*((1.70*sin(theta)+3.5*sqrt(alpha^2-0.25*sin(theta))+2*(alpha-sqrt(1-(p)^2))));
dy=A+B+C+D+E+F;
end
% end

Sign in to comment.

Tags

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!