Message Boards Message Boards

Non-linear transmission line: NDSolve error

Posted 2 years ago

Please I want to solve this equation in Mathematica

Clear["Global '*"];

(*nombre de section*)
NN = 11;

(*def du temps*)
tmax = 6.0;

(*def des parametres pricipaux de la ligne*)
fN = 1.0;                                                             \
                                (*facteur de normalisation*)
Rs = 50.0;                                                            \
                            (* resistance de la source en \
\[CapitalOmega]*)
Rescimpair = 
  Table[Rc[k] = 0.0, {k, 1, NN, 
    2}];           (*resistance du condensateur en \[CapitalOmega]*)
Rescpair = Table[Rc[k] = 0.0, {k, 2, NN, 2}];
Rl[NN] = 1100.0;                                                      \
                      (* resistance de la charge en \[CapitalOmega]*)


Reslimpair = 
  Table[Rl[k] = 0.0, {k, 1, NN, 
    2}];           (* resistance de l'inductance en \[CapitalOmega]*) \

Reslpair = Table[Rl[k] = 0.0, {k, 2, NN, 2}];

(*tension d'entrée*)
t0 = 0.0; ts = 1.7 + t0; tflat = 300.0; tfall = 1.7; tc = ts + tflat;
td = tc + tfall; a0 = -10.0*10^3; {ts, tflat, tfall, td - tc}
Vs[t_] := Which[t <= t0, 0, t > t0  && t <= ts, a0 *(t - t0)/(ts - t0),
  t > ts && t <= tc, a0, t > tc && t <= td, a0*((-t + td)/(td - tc)), 
  t > td, 0]
TabVs = Table[{t, Vs[t]/1000}, {t, 0, tmax, 0.01}];
VoltVs = ListPlot[TabVs, PlotRange -> All, 
  PlotStyle -> {AbsoluteThickness[1.4], RGBColor[0, 0, 1], 
    Thickness[0.008]}, FrameLabel -> {"Time(ns)", "Voltage(kV)"}, 
  Joined -> True, GridLines -> {Automatic, Automatic}, 
  FrameTicks -> {Automatic, Automatic}, Frame -> True, 
  LabelStyle -> Directive[Black, 17]]

(*transformation de fourier*)
np = 256; tfourier = 1000;
pulse = Table[Vs[t]/1000, {t, 0, tfourier, 0.05}];
datafin = Table[{f/(1.2), Abs[Fourier[pulse]][[f]]}, {f, 1, 60, 1}];
pVs = ListPlot[datafin, PlotRange -> {{0, 50}, All}, 
  PlotStyle -> {AbsoluteThickness[1.8], RGBColor[0, 1, 0], 
    Thickness[0.008]}, FrameLabel -> {"Frequence(Mhz)", "Voltage(V)"},
   Joined -> True, GridLines -> {Automatic, Automatic}, 
  FrameTicks -> {Automatic, Automatic}, Frame -> True, 
  LabelStyle -> Directive[Black, 17]]
NIntegrate[Vs[t]/1000, {t, 0, 200}]

(*non linear capacitor*)
m = 0.0; V0 = 0.7; Cs0 = 0.000095;
Cv[k_][t_] := which[v[k][t] > -V0, Cs0/p, v[k][t] <= -V0, 10.0]
Cv[k_][t_] := Cs0/(1.0 + 1*v[k][t]/V0)^m
Ca[Va_] := Cs0/(1.0 + 1*Va/V0)^m 
Cv[NN][t] := 10000000.0;
Plot[{Ca[Va]}, {Va, 0, 15}, AxesLabel -> {"Va", "Ca"}, 
 PlotRange -> {0, Cs0}, PlotStyle -> {Thickness[.01], Red}, 
 TicksStyle -> Directive["Black", 14], 
 AxesStyle -> {{Thick, Black}, {Thick, Black}}, 
 AxesLabel -> {Style["t", Black, Italic, 30], 
   Style["x,Vin", Black, Italic, 30]}, Frame -> False, 
 PlotRange -> All]

(*non linear inductance*)
L0 = 465.0; La = 4.65; Is = 3.76;
Ls[k_][t_] := (L0 - La)*(Sech[i[k][t]/Is]^2) + La
La1[Ia_] := (L0 - La)*(Sech[(Ia/Is)]^2) + La
Plot[{La1[Ia]}, {Ia, 0, 100}, AxesLabel -> {"Ia", "La1"}, 
 PlotRange -> {0, L0*1.2}, PlotStyle -> {Thickness[.01], Orange}, 
 TicksStyle -> Directive["Black", 14], 
 AxesStyle -> {{Thick, Black}, {Thick, Black}}, 
 AxesLabel -> {Style["t", Black, Italic, 30], 
   Style["x,Vin", Black, Italic, 30]}, Frame -> False, 
 PlotRange -> All]
(*Ls[NN][t]:=0.00000001*)
(*Ls[k_][t_]:=280.0*)

(*recursive fonction for current and charge*)
Denq = Table[Ls[k][t]*Cv[j][t], {j, 1, NN}];
Multi = Table[Rc[j]/Ls[k][t], {j, 1, NN}];
Ip[0][0] = 0.;
Qs[0][0] = 0.;
Qs[k_][t] := 
  Qs[k - 1][t] + 
   Sum[Cv[j][t]*v[j][t]/Denq[[k]], {j, k, NN}, {t, 0., tmax}];
Ip[k_][t] := 
  Ip[k - 1][t] + Multi[[k]]*Sum[i[j][t], {j, k, NN}, {t, 0., tmax}];

(* general equation for current and charge*)
eqs = Table[{i[k]'[t] == 
     Vs[t]/Ls[k][t] - (Rs/Ls[k][t])*
       Sum[i[j][t], {j, 1, NN}] - (Rl[k]/Ls[k][t] )*i[k][t] - Ip[k] - 
      Qs[k], v[k]'[t] == (i[k][t] - i[k + 1][t])/Cv[k][t]}, {k, 1, 
    NN}];
eqfinal = Flatten[eqs];

(*condition initiale*)
initial1 = 
  Flatten[Table[{i[k][0] == 0., v[k][0] == 0., Ip[k][0] == 0., 
     Qs[k][0] == 0.}, {k, 1, NN}]];
Vlist = Flatten[Table[{v[k][t], i[k][t]}, {k, 1, NN}]];
sol = NDSolve[Join[eqfinal, initial1], Vlist, {t, 0., tmax}, 
   MaxSteps -> Infinity];
sol1 = Flatten[sol];
inputiv = 
  Table[{i[k][t_] = i[k][t] //. sol1, 
    v[k][t_] = v[k][t] //. sol1}, {k, 1, NN}];
outiv = Flatten[inputiv];
V[0][t_] := Vs[t] - Rs*Sum[i[j][t], {j, 1, NN}];
V[NN][t_] := Rl[NN]*i[NN][t] + Ls[NN][t]*(i[NN][t] - i[NN + 1][t]);
Pload[t_] := V[NN][t] *i[NN][t]
Pint[t_] := Vs[t] *i[1][t]
TabVs1 = Table[{t, V[NN ][t]/1000}, {t, 0, tmax, 0.01}];
VoltFim = 
 ListPlot[TabVs1, PlotRange -> All, 
  PlotStyle -> {AbsoluteThickness[1.4], RGBColor[0, 1, 0], 
    Thickness[0.008]}, FrameLabel -> {"Time(ns)", "Voltage(kV)"}, 
  Joined -> True, GridLines -> {Automatic, Automatic}, 
  FrameTicks -> {Automatic, Automatic}, Frame -> True, 
  LabelStyle -> Directive[Black, 17]]
Show[VoltFim, VoltVs]

There are many problems with your code.

There is a lowercase which. Anyway, you should use Piecewise, and not Which, for symbolic math.

The equation

v[k]'[t] == (i[k][t] - i[k + 1][t])/Cv[k][t]

involves i[k+1]. When k ranges from 1 to NN the last equation involves i[NN+1].

Qs[0][t] is not defined, but is needed for Qs[1][t]

The definition

Ip[k_][t] := 
 Ip[k - 1][t] + Multi[[k]]*Sum[i[j][t], {j, k, NN}, {t, 0., tmax}]

uses t as both independent variable and iterator. What does it mean?

POSTED BY: Gianluca Gorni
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract