148 LETV(J)=m<J)-)-<ft<J)+F(J)*H>/$sHEXT J

158 GOSUB 488:F0R J=1 TO N!LETC<J)=F<J)iiiH

168 LETV<J)=u<J)+A<J)/8+.375*C(J)!HEXT J

178 LETX=X+H/6: GOSUB 488s FOR J=l TO Ы

188 LETB<J)=F<J)*H!LETV<J)=u<J)+A<J)/2-1.5iiiC<;j>+2iiiB<:j>

138 NEXT JS LETX=X+H-2! GOSUB 488

280 FOR J=l TO H:LETE<J)=F<J)*H

210 LETV<J)=u<J) + <:ft<J)+4iiiIi<J>+E<J))/6

220 LETE2=ABS<-2*ft<J)+3itiC<J)-8*J)<J>+E(J))/30

2.30 IF E2<=E1 THEN 250

240 LETE3=l!G0T0 266

250 IF E2<El/2e THEN LETD=B+1

260 NEXT JsIF E3=0 THEN 230

270 LETX=X-HsFOR J=l TO N!LETV<J)=u<J):NEXT J 280 LETH=H/2sG0T0 118 230 IF b=N THEN LETH=H+H

300 PRINT!F1.2PT=X,u BX=A7/u 3=v<1>,u C=V<2>

310 FOR J=l TO N!LETU<J)=V<J)!NEXT JsGOTO 118

400 IF X<=F THEN LETB7=u*<1-EXP<-X/T))

410 LETA7=B7!IF X>F THEN 1£Тй7=В7жЕХР<-<Х-Р>/Т>

420 LETIiS=A9+B9*EXP<-C9«V<2»!LETES=B9+E9*exp<:-F3*V<2)>

440 LETFS=mi< 1+HTN<S«<V< 1 )-u)/I))

458 LETV=F8*<1-exp<:-P*s*V(2)/F8))sLETA8=<:A7-V( 1)

470 LETBS= < < V< 3 ) -V> /-ES-ASz-C ) /- < 1 /-0+1 /B8+1 /-ES )

480 LETC8=V<3)-V-B8SLETF(1>=<A8+BS)/CSLETF(2)=C8/E8

490 LETF<:3)=<:E-V<3>*Z-V<2> >/As returns enb

Пример. Для !,= П A, Us= 12.6 В, S=2,72. А/В, K=QJ. С, =3,5-10-" Ф, C,2„„„ = 8.10-= Ф, Л = 3,5-10-" Ф, В = =0,088 1/В, С22м„„=1,3-10-" Ф, С = =6,59-10-" Ф, D=0,33 1/В (параметры транзистора КП913А), Ес = 45 В, Rc=5 Ом, £.= 1.10- Ф, t/„ex = 20 В, Твх = 3-10- с, <„=!4-10- с, йо=ЫО- с (начальный шаг), Q=l-10 = (погрешность вычислений) и to = 0 получаем следующие данные:

двух дифференциальных уравнений di Е~Ш - и du i-I {и)

где 1{U) - зависимость тока нелинейного двухполюсника от приложенного напряжения и. Отметим, что часть общей емкости С отражает инерционность двухполюсника.

Зависимости «вх (О. "зи (О и «СИ {t), рассчитанные по программе 7.61, показаны на рис. 7.22, причем на графике зависимости UcK{t) отмечены точки - результаты вычислений. Из расположения этих точек отчетливо виден характер автоматического изменения шага интегрирования дифференциальных уравнений - на крутых участках Uc„{f) шаг автоматически уменьшается, а на пологих - возрастает. В общей сложности для построения этих кривых понадобилось 72 цикла вычислений. Если бы использовался обычный метод Рунге - Кутта с шагом ft = 0,25-10- с (такой шаг нужен для получения погрешности Ql-10-=), то при общем временном интервале 30-10- с потребовалось бы 120 циклов вычислений. При увеличении протяженности пологих участков зависимости {Уси(<) выигрыш от применения метода Рунге - Кутта - Мерсона становится еще заметнее.

Моделирование нелинейной автоколебательной цепи (рис. 7.23) с двухполюсником, имеющим. Л-образную вольт-амперную характеристику, сводится к решению системы из

Допустим, что А/-образная вольт-амперная характеристика двухполюсника (в нашем случае туннельного диода ЗИ202К) задана графически (рис. 7.24). Аналитически ее можно представить полиномом 5-й степени, записанным по схеме Горнера

I (й) = ((((05И--04) и+аз) и + аг) и +

-fo,) и + ао. (7.21)

Для нахождения коэффициентов оо-as полинома составим таблицу значений и и Г.

Используя любую из описанных ранее программ для полиномиальной интерполяции и аппроксимации при п = 5 и ft=At/=0,2 В, находим: ао = О, d =746,333, 02 =-3420,833, 03 = 5734,375. 04=-4166,667 и 05=1119,792. Кривая, построенная по формуле (7.21), практически точно совпадает с исходным

.....

30 20

Рис. 7.22. Расчетные зависимости "sm(0 и «„(О ялн ключа на мощном f-МДП-транзисторе,

полученные при его моделировании

Рнс. 7.23. Нелинейная автоколебательная схема с двухполюсником, имеющим /У-образную вычислительную характеристику

/-1-Ш)

Предельный цикл I колебаний

-ZOV

Рис. 7.24. Фазовый портрет автоколебаний, построенный Б плоскости вольт-амперной характеристики

графиком зависимости /(ы), что свидетельствует о высокой точности аппроксимации.

Для решения приведенной выше системы воспользуемся методом Рунге - Кутта -

Мерсона, дополнительно организовав ввод исходны! данных и подпрограмму для вычисления / (и), di/dt = FO и du/dt = F\ при 1 = Л0 и и = А\.

10 PR i NTМОДЕЛИРОВАНИЕ РЕЛАКСАЦИОННОГО ГЕНЕРАТОРА

20 printha ПРИБОРЕ С н-ОБРАЗНОй В.А.Х.

25 letn=2:bim VaO/w*:N>;A*:N),c<h>TlKN>.e<h>TFCN>

30 printВВЕДИТЕ КОЭФФИЦИЕНТЫ АППРОКСИМАЦИИ Е.а.Х.

40 input а0/ AIjА2/ a3/ а4,а5q0, q1, 02/ 03/ 04/ 65

50 inputВВЕДИТЕ ПАРАМЕТРЫ СХЕМЫ e/r/c/LE,Z/c/L

60 inputЗАДАЙТЕ НАЧАЛЬНЫЙ ШАГ н0=н

70 inputЗАДАЙТЕ ПОГРЕШНОСТЬ EPS=e1

90 let vo > =0 5 letu а > =0 S letv < 2 > =0 i LETIK 2 > =0: letx=e

110 lete3=0:6osub 400!letd=0!for J=l ТОН

120 leta<;j>=f<.J>3iiH!letv<;.j>=lAKJ>+ft<J>--3!hext J

130 letx=x+h.3!60sue 400!for J=l to h

140 letv с j > =lj.i с j > + < A < j > +f < j > ЖН > -6! next j

150 60sub 400!for J=l to N!letccj>=fcj>*h

160 letv<;j>=u<;j>+ft*:.j>8+.375.*c<;j>!next j

170 letx=x+h.6s60SUB 400! for J=l to N

180 LETIi<.J>=F<j>*h!lety<J>=IKj>+a<j>/2-l.5.*c<;J>+2*IK.J>

190 next j!letx=x+h/2!60sub 488

200 for j=l to N!lete<;j>=f(j)*h

210 letvcj>=w*:j>+<:a<:.j>+4*d<;j>+e<;j>>.6

220 lete2=abs ( -2.*a < j >+9.*c с j > -8«d с j ) +е с j > > -38

230 if e2<=E1 then 258

240 lete3=l:60t0 268

250 if E2<e 1-20 then letij=Ii+l

260 next j!if E3=0 then 296

270 letx=x-h!for j=l to N!letv<:.j>=i.Kj)!hext j 280 leth=h/2!60t0 118 290 if d-N then leth=h+h

300 print!f1.2!t=x/u=v<:2>/ir=Va>*lE3

310 letw<1>=v<1>!letw<:2>=V<2):60t0 118

400 letv=V*:2>!leti=<; ii cq5*V+e4>.*v+Q3>.*v+G2>«v+Ul>*v+u0

410 LETFa>=<:e-va>*2-v<:2>>/L

420 letf с 2 >=с V < 1 > -1 1e3 > -c !/return: end

£ t,HD

Рис. 7.25. Расчетные временные завнсимости u(l) и /(/) для автоколебательной схемы на туннельном

диоде

Пример. Для указанных значений Со-о,. £=0,5 В, /?=5 Ом, С= 100-10-= Ф, Z, = 50-10- Гн, йо=!-!0-° с (начальный шаг). Q = l-10-= (погрешнсх;ть вычислений) и <о=0 будем получать следующие результаты:

На рис. 7.24 показана фазовая траектория автоколебаний при нулевых начальных условиях, а на рис. 7.25 - форма импульсов напряжения u{f) и тока i{t\. Аналогичные расчетным данные получены при экспериментальном исследовании данной схемы.