All I think is to get a sine wave function when plotting T vs A(1), but I didn't get any output. help me to modify this program

조회 수: 1 (최근 30일)
format long
a = 0;
b = 4;
h = 0.1
V =1E-5
I1 = 0.2;
I2 = 0.8;
o = 3E5;
tc = 30E-9;
tf = 230E-6;
a1 = 0.1;
a2 =0.1;
P1 = 0.2;
P2 =0.2;
k= 0.17;
A1(1) = 0;
A2(1) = 0;
G1(1) = ;
G2(1) = 0;
O(1) = 0;
n = (b-a)/h ;
t = a + (0:n)*h;
func1 = @(G1,A1,A2,O) (1/tc).*(G1- a1).*A1 +(k./tc).*A2.*cos(O);
func3 = @(G2,A2,A1,O) (1/tc).*(G2- a2).*A2 +(k./tc).*A1.*cos(O);
func2 = @(G1) (1/tf).*(P1 - G1.*(I1+1));
func4 = @(G2) (1/tf).*(P2 - G2.*(I2+1));
func5 = @(A1,A2,O) o - (k./tc).*((A1./A2) + (A2./A1)).*sin(O);
for i = 1:n
k1 = func1(G1(i),A1(i),A2(i),O(i));
l1 = feval(func2, G1(i));
m1 = feval(func3, G2(i),A1(i),A2(i),O(i));
n1 = feval(func4, G2(i));
q1 = feval(func5,A1(i),A2(i),O(i));
k2 = feval(func1, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h, G1(i)+(l1/2)*h, O(i)+(q1/2)*h);
l2 = feval(func2, G1(i)+(l1*h/2));
m2 = feval(func3, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h, G2(i)+(n1/2)*h, O(i)+(q1/2)*h);
n2 = feval(func4, G2(i)+(n1/2)*h);
q2 = feval(func5, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h, O(i)+(q1/2)*h );
k3 = feval(func1, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h, G1(i)+(l2/2)*h, O(i)+(q2/2)*h);
l3 = feval(func2, G1(i)+(l2*h/2));
m3 = feval(func3, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h, G2(i)+(n2/2)*h, O(i)+(q2/2)*h);
n3 = feval(func4, G2(i)+(n2/2)*h);
q3 = feval(func5, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h, O(i)+(q2/2)*h );
k4 = feval(func1, A2(i)+(m3*h), A1(i)+(k3*h), G1(i)+(l3*h), O(i)+(q3*h));
l4 = feval(func2, G1(i)+(l3*h));
m4 = feval(func3, A2(i)+(m3*h), A1(i)+(k3*h), G2(i)+(n3*h), O(i)+(q3*h));
n4 = feval(func4, G2(i)+(n3*h));
q4 = feval(func5, A2(i)+(m3*h), A1(i)+(k3*h), O(i)+(q3*h) );
A1(i+1) = A1(i) + (h/6)*(k1+(2*(k2+k3))+k4);
G1(i+1) = G1(i) + (h/6)*(l1+(2*(l2+l3))+l4);
A2(i+1) = A2(i) + (h/6)*(m1+(2*(m2+m3))+m4);
G2(i+1) = G2(i) + (h/6)*(n1+(2*(n2+n3))+n4);
O(i+1) = O(i) + (h/6)*(q1+(2*(q2+q3))+q4);
end
subplot 321
plot(t,A1)
subplot 322
plot(t,A2)
subplot 323
plot(t,G1)
subplot 324
plot(t,G1)
subplot 325
plot(t,O)
  댓글 수: 3
SAHIL SAHOO
SAHIL SAHOO 2022년 7월 3일
hey sam,
thanks for reply,
and this program I suppose to solve the 5 couple oscillator and the ouput should be sinusoidal wave between A(1) vs t.
these are equations that have to be solve by the ODE method.
dpb
dpb 2022년 7월 3일
Have you verified the calculation of your functions at the origin?
>> i = 1
k1 = func1(G1(i),A1(i),A2(i),O(i))
l1 = feval(func2, G1(i))
m1 = feval(func3, G2(i),A1(i),A2(i),O(i))
n1 = feval(func4, G2(i))
q1 = feval(func5,A1(i),A2(i),O(i))
i =
1
k1 =
0
l1 =
869.5652
m1 =
0
n1 =
869.5652
q1 =
NaN
>>
shows your func5 returns NaN initially and since plot() doesn't show NaN values, only the origin 0 point will be plotted. Similar problems exist for the other functions as well and your integration grows without bounds for those two for which it doesn't return NaN.
Need to use the debugger and step through and see where your formulation went wrong...
BTW, you can simplify the coding -- there's no need to use feval here;, just use the anonymous function handles with the argument lists as you did in the very first call for k1 everywhere--not sure why you would have changed after that line for l1, m1, ...
k1 = func1(G1(i),A1(i),A2(i),O(i));
l1 = func2(G1(i));
m1 = func3(G2(i),A1(i),A2(i),O(i));
n1 = func4(G2(i));
q1 = func5(A1(i),A2(i),O(i));
...
Unless this is homework and required to write the integration explicitly as part of the assignment, would be easier to use the built-in ODE solvers in MATLAB.

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

채택된 답변

VBBV
VBBV 2022년 7월 3일
format short
a = 0;
b = 4;
h = 0.1;
V =1E-5;
I1 = 0.2;
I2 = 0.8;
o = 3.5; % what is this variable ?
tc = 3.0; %
tf = 2.3; %
a1 = 0.1;
a2 =0.1;
P1 = 0.2;
P2 =0.2;
k= 0.17;
% initial conditions
A1(1) = 0.7;
A2(1) = 0.5;
G1(1) = 0.5;
G2(1) = 0.5;
O(1) = 10;
n = (b-a)/h ;
t = a + (0:n)*h;
func1 = @(G1,A1,A2,O) (1/tc).*(G1- a1).*A1 +(k./tc).*A2.*cos(O);
func3 = @(G2,A2,A1,O) (1/tc).*(G2- a2).*A2 +(k./tc).*A1.*cos(O);
func2 = @(G1) (1/tf).*(P1 - G1.*(I1+1));
func4 = @(G2) (1/tf).*(P2 - G2.*(I2+1));
func5 = @(A1,A2,O) O - (k./tc).*((A1./A2) + (A2./A1)).*sin(O);
for i = 1:n
k1 = feval(func1,G1(i),A1(i),A2(i),O(i));
l1 = feval(func2, G1(i));
m1 = feval(func3, G2(i),A1(i),A2(i),O(i));
n1 = feval(func4, G2(i));
q1 = feval(func5,A1(i),A2(i),O(i));
k2 = feval(func1, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h, G1(i)+(l1/2)*h, O(i)+(q1/2)*h);
l2 = feval(func2, G1(i)+(l1*h/2));
m2 = feval(func3, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h, G2(i)+(n1/2)*h, O(i)+(q1/2)*h);
n2 = feval(func4, G2(i)+(n1/2)*h);
q2 = feval(func5, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h, O(i)+(q1/2)*h);
k3 = feval(func1, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h, G1(i)+(l2/2)*h, O(i)+(q2/2)*h);
l3 = feval(func2, G1(i)+(l2*h/2));
m3 = feval(func3, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h, G2(i)+(n2/2)*h, O(i)+(q2/2)*h);
n3 = feval(func4, G2(i)+(n2/2)*h);
q3 = feval(func5, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h, O(i)+(q2/2)*h );
k4 = feval(func1, A2(i)+(m3*h), A1(i)+(k3*h), G1(i)+(l3*h), O(i)+(q3*h));
l4 = feval(func2, G1(i)+(l3*h));
m4 = feval(func3, A2(i)+(m3*h), A1(i)+(k3*h), G2(i)+(n3*h), O(i)+(q3*h));
n4 = feval(func4, G2(i)+(n3*h));
q4 = feval(func5, A2(i)+(m3*h), A1(i)+(k3*h), O(i)+(q3*h));
A1(i+1) = A1(i) + (h/6)*(k1+(2*(k2+k3))+k4);
G1(i+1) = G1(i) + (h/6)*(l1+(2*(l2+l3))+l4);
A2(i+1) = A2(i) + (h/6)*(m1+(2*(m2+m3))+m4);
G2(i+1) = G2(i) + (h/6)*(n1+(2*(n2+n3))+n4);
O(i+1) = O(i) + (h/6)*(q1+(2*(q2+q3))+q4);
end
subplot 321
plot(t,A1)
subplot 322
plot(t,A2)
subplot 323
plot(t,G1)
subplot 324
plot(t,G1)
subplot 325
plot(t,O)
Set intial conditions in a suitable to get desired result.
  댓글 수: 2
VBBV
VBBV 2022년 7월 3일
Check also using ode45 or similar builtin functions and varying input parameter values, for comparison purposes.
SAHIL SAHOO
SAHIL SAHOO 2022년 7월 5일
"Unable to perform assignment because the left and right sides have a different number of elements.
Error in laser2 (line 54)
A1(i+1) = A1(i) + (h/6)*(k1+(2*(k2+k3))+k4); "
I'm having this problem when I change
func2 = @(G1) (1/tf).*(P1 - G1.*(A1.*A1+1));
func4 = @(G2) (1/tf).*(P2 - G2.*(A2.*A2+1));
why this is happening, let me know how to solve this one?

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

추가 답변 (0개)

카테고리

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

Community Treasure Hunt

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

Start Hunting!

Translated by