while solving a couple ordinary differential equation, I'm getting an error of Unable to perform assignment because the left and right sides have a different number of element

4 views (last 30 days)
output:-
Unable to perform assignment because the left and right sides have a different number of element.
Error in laser2 (line 52)
A1(i+1) = A1(i) + (h/6)*(k1+(2*(k2+k3))+k4);
input:-
format long
a = 0;
b = 30;
h = 0.1;
o = 3.5;
tc = 3.0;
tf = 2.3;
a1 = 0.1;
a2 = 0.1;
P1 = 0.2;
P2 = 0.2;
k= 0.9;
V = 1;
A1(1) = 10E-4;
A2(1) = 10E-5;
G1(1) = 10E-6;
G2(1) = 10E-6;
O(1) = 0;
n = (b-a)/h ;
t = ((a) + (0:n)*h).*10E-6;
func1 = @(G1,A1,A2,O) (1/tc).*(G1- a1).*(A1./V) +(k./tc).*(A2./V).*cos(O);
func3 = @(G2,A2,A1,O) (1/tc).*(G2- a2).*(A2./V) +(k./tc).*(A1./V).*cos(O);
func2 = @(G1) (1/tf).*(P1 - G1.*(A1.^2+1));
func4 = @(G2) (1/tf).*(P2 - G2.*(A2.^2+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 211
plot(t,A1)
grid on
subplot 212
plot(t,O)
it maybe due to the func2 = @(G1) (1/tf).*(P1 - G1.*(A1.^2+1));
func4 = @(G2) (1/tf).*(P2 - G2.*(A2.^2+1));
this two code, but I couldn't understand what's the problem in it.

Accepted Answer

Johan
Johan on 6 Jul 2022
func2 and func4 uses A1 and A2 respectively which means they are using the whole arrays, either call A1 and A2 in the func definition (see example below) or set it to A1(1) and A2(1) if you only want to use the first term.
On a side note it's considered bad practice to add 1 2 3 .... to your variable names. It makes your code hard to understand and debug there is a post that explain why here : https://fr.mathworks.com/support/search.html/answers/304528-tutorial-why-variables-should-not-be-named-dynamically-eval.html
%change some variable names so that they make sense
time_start = 0;
time_end = 30;
time_step = 0.1;
o = 3.5;
tc = 3.0;
tf = 2.3;
a1 = 0.1;
a2 = 0.1;
P1 = 0.2;
P2 = 0.2;
k= 0.9;
V = 1;
n_step = (time_end-time_start)/time_step ;
time = ((time_start) + (0:n_step)*time_step).*10E-6;
func1 = @(G1,A1,A2,O) (1/tc).*(G1- a1).*(A1./V) +(k./tc).*(A2./V).*cos(O);
func3 = @(G2,A2,A1,O) (1/tc).*(G2- a2).*(A2./V) +(k./tc).*(A1./V).*cos(O);
func2 = @(G1,A1) (1/tf).*(P1 - G1.*(A1.^2+1));
func4 = @(G2,A2) (1/tf).*(P2 - G2.*(A2.^2+1));
func5 = @(A1,A2,O) o - (k./tc).*((A1./A2) + (A2./A1)).*sin(O);
%Initialize arrays
A1 = zeros(n_step,1);
A2 = A1;
G1 = A1;
G2 = A1;
O = A1;
A1(1) = 10E-4;
A2(1) = 10E-5;
G1(1) = 10E-6;
G2(1) = 10E-6;
O(1) = 0;
for i_loop = 1:n_step
k1 = func1(G1(i_loop),A1(i_loop),A2(i_loop),O(i_loop));
l1 = func2( G1(i_loop),A1(i_loop));
m1 = func3( G2(i_loop),A1(i_loop),A2(i_loop),O(i_loop));
n1 = func4( G2(i_loop),A2(i_loop));
q1 = func5(A1(i_loop),A2(i_loop),O(i_loop));
k2 = func1( A2(i_loop)+(m1/2)*time_step, A1(i_loop)+(k1/2)*time_step, G1(i_loop)+(l1/2)*time_step, O(i_loop)+(q1/2)*time_step);
l2 = func2( G1(i_loop)+(l1*time_step/2), A1(i_loop)+(k1/2)*time_step);
m2 = func3( A2(i_loop)+(m1/2)*time_step, A1(i_loop)+(k1/2)*time_step, G2(i_loop)+(n1/2)*time_step, O(i_loop)+(q1/2)*time_step);
n2 = func4( G2(i_loop)+(n1/2)*time_step, A2(i_loop)+(m1/2)*time_step);
q2 = func5( A2(i_loop)+(m1/2)*time_step, A1(i_loop)+(k1/2)*time_step, O(i_loop)+(q1/2)*time_step);
k3 = func1( A2(i_loop)+(m2/2)*time_step, A1(i_loop)+(k2/2)*time_step, G1(i_loop)+(l2/2)*time_step, O(i_loop)+(q2/2)*time_step);
l3 = func2( G1(i_loop)+(l2*time_step/2), A1(i_loop)+(k2/2)*time_step);
m3 = func3( A2(i_loop)+(m2/2)*time_step, A1(i_loop)+(k2/2)*time_step, G2(i_loop)+(n2/2)*time_step, O(i_loop)+(q2/2)*time_step);
n3 = func4( G2(i_loop)+(n2/2)*time_step, A2(i_loop)+(m2/2)*time_step);
q3 = func5( A2(i_loop)+(m2/2)*time_step, A1(i_loop)+(k2/2)*time_step, O(i_loop)+(q2/2)*time_step );
k4 = func1( A2(i_loop)+(m3*time_step), A1(i_loop)+(k3*time_step), G1(i_loop)+(l3*time_step), O(i_loop)+(q3*time_step));
l4 = func2( G1(i_loop)+(l3*time_step), A1(i_loop)+(k3*time_step));
m4 = func3( A2(i_loop)+(m3*time_step), A1(i_loop)+(k3*time_step), G2(i_loop)+(n3*time_step), O(i_loop)+(q3*time_step));
n4 = func4( G2(i_loop)+(n3*time_step),A2(i_loop)+(m3*time_step));
q4 = func5( A2(i_loop)+(m3*time_step), A1(i_loop)+(k3*time_step), O(i_loop)+(q3*time_step));
A1(i_loop+1) = A1(i_loop) + (time_step/6)*(k1+(2*(k2+k3))+k4);
G1(i_loop+1) = G1(i_loop) + (time_step/6)*(l1+(2*(l2+l3))+l4);
A2(i_loop+1) = A2(i_loop) + (time_step/6)*(m1+(2*(m2+m3))+m4);
G2(i_loop+1) = G2(i_loop) + (time_step/6)*(n1+(2*(n2+n3))+n4);
O(i_loop+1) = O(i_loop) + (time_step/6)*(q1+(2*(q2+q3))+q4);
end
subplot 211
plot(time,A1)
grid on
subplot 212
plot(time,O)

More Answers (2)

KSSV
KSSV on 6 Jul 2022
This function:
func2 = @(G1) (1/tf).*(P1 - G1.*(A1.^2+1));
should be:
func2 = @(G1,A1) (1/tf).*(P1 - G1.*(A1.^2+1));

Walter Roberson
Walter Roberson on 6 Jul 2022
The posted code executes for me
format long
a = 0;
b = 30;
h = 0.1;
o = 3.5;
tc = 3.0;
tf = 2.3;
a1 = 0.1;
a2 = 0.1;
P1 = 0.2;
P2 = 0.2;
k= 0.9;
V = 1;
A1(1) = 10E-4;
A2(1) = 10E-5;
G1(1) = 10E-6;
G2(1) = 10E-6;
O(1) = 0;
n = (b-a)/h ;
t = ((a) + (0:n)*h).*10E-6;
func1 = @(G1,A1,A2,O) (1/tc).*(G1- a1).*(A1./V) +(k./tc).*(A2./V).*cos(O);
func3 = @(G2,A2,A1,O) (1/tc).*(G2- a2).*(A2./V) +(k./tc).*(A1./V).*cos(O);
func2 = @(G1) (1/tf).*(P1 - G1.*(A1.^2+1));
func4 = @(G2) (1/tf).*(P2 - G2.*(A2.^2+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 211
plot(t,A1)
grid on
subplot 212
plot(t,O)

Community Treasure Hunt

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

Start Hunting!

Translated by