while solving coupled differential equation getting this problem.

"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 in this program

format short

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= 1.2;
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)

why this is happening, let me know how to solve this one?

[Torsten]

Works with MATLAB's R2022a.

[SAHIL SAHOO]

why this is terminate in 3*10^-4, could you explain this? and how to resolve this and go further ?

because by chnaging b= 40 still it's terminate in 3

[SAHIL SAHOO]

it's still showing error after I ran the code in matlab R2022a

Answer 1

[Ishan]

The dimension mismatch you get is because A1(i+1) = A1(i) + (h/6)*(k1+(2*(k2+k3))+k4); in this equation on the LHS you have allocated A1 to be of size 1 times 1 while the RHS is 1 times 301(that too if the loop iterates all the way through without throwing an error).

So basically, you assigned the value of A1 as 10e-4 which is a ‘double’ value of 1*1 size. Now when you index within the ‘for’ loop what is happening is that for each iteration a new Value is being appended to all the variables (thereby increasing their dimension on each iteration). As far as the graph stopping before the expected end point, that is because your loop runs for 300 iterations so the value of the variables you want to plot are calculated for the corresponding number of times. A simple fix would be to plot for over 300 values only. Although the exact coupled differential equation is unknown to me to be able to validate the results, this code should be able to resolve your current issue: -

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= 1.2;
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,A1,A2,O);
    %do not use l1(i), instead you can store value of l1 in each iteration
    %in that same variable
    l1 = feval(func2, G1);
    m1 = feval(func3, G2,A1,A2,O);
    n1 = feval(func4, G2);
    q1 = feval(func5,A1,A2,O);
    
    k2 = feval(func1, A2+(m1/2)*h, A1+(k1/2)*h,  G1+(l1/2)*h,  O+(q1/2)*h);
    l2 = feval(func2, G1+(l1*h/2));
    m2 = feval(func3, A2+(m1/2)*h, A1+(k1/2)*h,  G2+(n1/2)*h,  O+(q1/2)*h);
    n2 = feval(func4, G2+(n1/2)*h);
    q2 = feval(func5, A2+(m1/2)*h, A1+(k1/2)*h, O+(q1/2)*h);
    
    k3 = feval(func1, A2+(m2/2)*h, A1+(k2/2)*h,  G1+(l2/2)*h,  O+(q2/2)*h);
    l3 = feval(func2, G1+(l2*h/2));
    m3 = feval(func3, A2+(m2/2)*h, A1+(k2/2)*h,  G2+(n2/2)*h,  O+(q2/2)*h);
    n3 = feval(func4, G2+(n2/2)*h);
    q3 = feval(func5, A2+(m2/2)*h, A1+(k2/2)*h, O+(q2/2)*h );
    
    k4 = feval(func1, A2+(m3*h), A1+(k3*h),  G1+(l3*h),  O+(q3*h));
    l4 = feval(func2, G1+(l3*h));
    m4 = feval(func3, A2+(m3*h), A1+(k3*h),  G2+(n3*h),  O+(q3*h));
    n4 = feval(func4, G2+(n3*h));
    q4 = feval(func5, A2+(m3*h), A1+(k3*h), O+(q3*h));
    
    size((h/6)*(k1+(2*(k2+k3))+k4));
    A1 = A1 + (h/6)*(k1+(2*(k2+k3))+k4);
    G1 = G1 + (h/6)*(l1+(2*(l2+l3))+l4);
    A2 = A2 + (h/6)*(m1+(2*(m2+m3))+m4); 
    G2 = G2 + (h/6)*(n1+(2*(n2+n3))+n4); 
    O  = O + (h/6)*(q1+(2*(q2+q3))+q4);
    %store A1 and O in separate variables and plot it w.r.t t (300 length)
    y(i) = O;
    x(i) = A1;
end
subplot 211
t = 1:1:300;
plot(t,x)
grid on
subplot 212
plot(t,y)

Note: You can also scale the axis as per your requirements