how to use for loop in equation

Tstart = 0.0;

Tend = 1000.0;

Nt = 20000;

dT = (Tend-Tstart)/Nt;

X0 = 0.0;

Y0 = 1.0;

Z0 = 0.0;

N = 20;

for A=0.1:1

B=0.2;

C=5.7;

T = zeros(Nt+1,1);

X = zeros(Nt+1,1);

Y = zeros(Nt+1,1);

Z = zeros(Nt+1,1);

a = zeros(N+1,1);

b = zeros(N+1,1);

c = zeros(N+1,1);

T(1) = 0.0;

X(1) = X0;

Y(1) = Y0;

Z(1) = Z0;

for j = 2:Nt+1

a(1) = X(j-1);

b(1) = Y(j-1);

c(1) = Z(j-1);

for k = 1:N

SC = 0.0;

for i= 0:k-1

SC = SC+ a(i+1)*c(k-i);

end

a(k+1) = (-b(k) - c(k))/k;

b(k+1) = (a(k)+A*b(k))/k ;

c(k+1) = (SC+B-C*(c(k)))/k ;

end

x = a(1);

y = b(1);

z = c(1);

for k = 2:N+1

x = x + a(k)*dT^(k-1);

y = y + b(k)*dT^(k-1);

z = z + c(k)*dT^(k-1);

end

end

T(j) = T(j-1) + dT;

X(j) = x;

Y(j) = y;

Z(j) = z;

end

plot(A,X)

how to plot the A vs X

[shiv gaur]

what is the problem in looping I am not able to plot

Answer 1

[Riccardo Scorretti]

Hi Shiv,

up to my understanding you want to execute a kind of parametric study, basing on a parameter A (and perhaps in future on B and C). If so, the point is that:

1. indeed, A iterates only one time with value A = 0.1,
2. most importantly, you was trying to plot X (= a vector) versus A (= a scalar). If the purpose is to plot the result of computation for several values of A, you should plot X as a function of T inside the outer loop, and hold on the graphic.

Perhaps, you meant to do this?

Tstart = 0.0;
Tend = 1000.0;
Nt = 20000;
dT = (Tend-Tstart)/Nt;
X0 = 0.0;
Y0 = 1.0;
Z0 = 0.0;
N = 20;

figure;                         % ***
for A=0.1:1
    B=0.2;
    C=5.7;
    
    T = zeros(Nt+1,1);
    X = zeros(Nt+1,1);
    Y = zeros(Nt+1,1);
    Z = zeros(Nt+1,1);
    a = zeros(N+1,1);
    b = zeros(N+1,1);
    c = zeros(N+1,1);
    
    T(1) = 0.0;
    X(1) = X0;
    Y(1) = Y0;
    Z(1) = Z0;
    
    
    for j = 2:Nt+1
        
        a(1) = X(j-1);
        b(1) = Y(j-1);
        c(1) = Z(j-1);
        for k = 1:N
        
            SC = 0.0;
            for i= 0:k-1
                SC = SC+ a(i+1)*c(k-i);
                
            end
            a(k+1) = (-b(k) - c(k))/k;
            b(k+1) = (a(k)+A*b(k))/k ;
            c(k+1) = (SC+B-C*(c(k)))/k ;
        end 
        
        
        
        x = a(1);
        y = b(1);
        z = c(1);
        for k = 2:N+1
            x = x + a(k)*dT^(k-1);
            y = y + b(k)*dT^(k-1);
            z = z + c(k)*dT^(k-1);
        end
    end
        
        
    T(j) = T(j-1) + dT;
    X(j) = x;
    Y(j) = y;
    Z(j) = z;
    
    plot(T,X);              % ***
    hold on;                % ***
end

% plot(A,X)