Your Y matrix is still wrong. And what about R? And, as @zdenek_hurak said, fix the slash.
You must look carefully over every line of your code, please. Every last thing must be correctly translated.
1 Like
I checked my code again and again . Thanks for the reply and reminding me about my silly mistake but the figures of MATLAB and Julia are not the same
using Plots
plotly()
Plots.PlotlyBackend()
currentcap=zeros(1,250001);
voltagecap=zeros(1,250001);
currentind=zeros(1,250001);
V= zeros(9,250004);
t_end = 250e-3;
dt = 1e-6;
t = 0:dt:t_end;
vp=20;
freq=50;
vs = vp*sin.(2*π*freq.* t')
L=1e-3;
R=10000;
C=1e-6;
Roff=10^10;
a=0; b=1; f=0; g=1; c=0; e=1; h=0; n=1;
for i=2:250001
Y = [dt/L -dt/L 0 0 1 0 0 0 0
-dt/L dt/L+1/Roff -c/Roff -a/Roff 0 1 0 0 -1
0 -c/Roff (1/R)+(C/dt)+1/Roff -1/R-C/dt 0 -1 -1 0 0
0 -a/Roff -1/R-C/dt 1/R+C/dt+1/Roff 0 0 0 1 1
1 0 0 0 0 0 0 0 0
0 a -a 0 0 b 0 0 0
0 0 -c 0 0 0 e 0 0
0 0 0 -f 0 0 0 g 0
0 -h 0 h 0 0 0 0 n];
x1=V[5,i]-currentind[1,i-1]; x2=currentind[1,i-1]-V[6,i]+V[9,i]; x3=V[6,i]-currentcap[1,i-1]+V[7,i]; x4=currentcap[1,i]-V[8,i]-V[9,i] ;
A =[x1
x2
x3
x4
vp*sin.(2*π*freq.*t'[1,i])
0
0
0
0];
x=A\Y;
V[:,i]=x;
currentcap[1,i]=(C/dt)*(V[3,i]-V[3,i-1]-V[4,i]+V[4,i-1]);
voltagecap[1,i]=V[3,i]-V[4,i];
currentind[1,i]=(dt/L)*(V[1,i]-V[2,i]+V[1,i-1]-V[2,i-1])+currentind[1,i-1];
if V[2,i]-V[3,i]>0 || V[6,i]>0
a=1;b=0;
elseif V[6,i]<=0
a=0; b=1;
else
a=0; b=1;
end
if -V[3,i]>0 || V[7,i]>0
c=1; e=0;
elseif V[7,i]<=0
c=0; e=1
else
c=0; e=1;
end
if V[4,i]>0 || V[8,i]>0
f=1; g=0;
elseif V[8,i]<=0
f=0; g=1
else
f=0; g=1;
end
if V[4,i]-V[2,i]>0 || V[9,i]>0
h=1; n=0;
elseif V[9,i]<=0
h=0; n=1
else
h=0; n=1;
end
end
voltagecap=voltagecap';
plot(t,voltagecap)
the plot is as follows:
but the figure of Matlab is like this:
Thanks for the reply,
I corrected the Matrix and the value of R but still the output of my code in Julia is not correct
Your linear equation inside the code is still not set up correctly. In your original Matlab code, you form a matrix Y and a column vector I. You then go ahead and solve the equation Yx=I for a vector x. In your code you do it by computing the inverse of the matrix as in x=inv(Y)*I, which is not a recommendable way from a numerical point of view, and therefore you were advised to use backslash. The same recommendation to avoid explicit computation of the inverse holds in Julia. And coincidently, the same syntax with backslash is used in Julia. You were even given the particular code: x = Y\I. Why did not you just use it?
Instead, in your new code I can observe that you just renamed your vector I to A and now you try to solve the set of linear equations by writing x=A\Y, which sends you back to your previous (incorrect) solution x = I\Y.
Perhaps there is some misunderstanding as for what the backslash means here. If you do not find the usage of the symbol intuitive, have a look at the documentation for backslash. If you want to solve Ax=b for x, given some matrix A and (column) vector b, you just write x=A\b. But this is identical to Matlab syntax.
2 Likes
I edited the Y\A but I don’t know why the code never enters the loop. It’s so strange. This is my first problem.
using Plots
plotly()
Plots.PlotlyBackend()
currentcap=zeros(1,250001);
voltagecap=zeros(1,250001);
currentind=zeros(1,250001);
V= zeros(9,250004);
t_end = 250e-3;
dt = 1e-6;
t = 0:dt:t_end;
vp=20;
freq=50;
vs = vp*sin.(2*π*freq.* t')
L=1e-3;
R=10000;
C=1e-6;
Roff=10^10;
c1=0;
a=0; b=1; f=0; g=1; c=0; e=1; h=0; n=1;
i=2;
for i=2:250001
c1=c1+1;
Y = [dt/L -dt/L 0 0 1 0 0 0 0
-dt/L dt/L+1/Roff -c/Roff -a/Roff 0 1 0 0 -1
0 -c/Roff (1/R)+(C/dt)+1/Roff -1/R-C/dt 0 -1 -1 0 0
0 -a/Roff -1/R-C/dt 1/R+C/dt+1/Roff 0 0 0 1 1
1 0 0 0 0 0 0 0 0
0 a -a 0 0 b 0 0 0
0 0 -c 0 0 0 e 0 0
0 0 0 -f 0 0 0 g 0
0 -h 0 h 0 0 0 0 n];
x1=V[5,i]-currentind[1,i-1]; x2=currentind[1,i-1]-V[6,i]+V[9,i]; x3=V[6,i]-currentcap[1,i-1]+V[7,i]; x4=currentcap[1,i]-V[8,i]-V[9,i] ;
A =[x1
x2
x3
x4
vp*sin.(2*π*freq.*t'[1,i])
0
0
0
0];
x=Y\A;
V[:,i]=x;
currentcap[1,i]=(C/dt)*(V[3,i]-V[3,i-1]-V[4,i]+V[4,i-1]);
voltagecap[1,i]=V[3,i]-V[4,i];
currentind[1,i]=(dt/L)*(V[1,i]-V[2,i]+V[1,i-1]-V[2,i-1])+currentind[1,i-1];
if V[2,i]-V[3,i]>0 || V[6,i]>0
a=1;b=0;
elseif V[6,i]<=0
a=0; b=1;
else
a=0; b=1;
end
if -V[3,i]>0 || V[7,i]>0
c=1; e=0;
elseif V[7,i]<=0
c=0; e=1
else
c=0; e=1;
end
if V[4,i]>0 || V[8,i]>0
f=1; g=0;
elseif V[8,i]<=0
f=0; g=1
else
f=0; g=1;
end
if V[4,i]-V[2,i]>0 || V[9,i]>0
h=1; n=0;
elseif V[9,i]<=0
h=0; n=1
else
h=0; n=1;
end
end
voltagecap=voltagecap';
plot(t,voltagecap)