nonlinear_dynamic_example_fdtd.m

non linear simulation techniques for circuits FDTD

% *************************************************************************
% This file allows to simulate the nonlinear dynamic circuit example
% described by the following equation: 
%
%                       Gvo(t)+dqnl(vo(t))/dt+inl(vo(t))=is(t)
%
% The method used for that will be FDTD (Finite Diferences in time domain)
% *************************************************************************

clear all;
close all;
clc;

%Constants definition
raphson_error=0.000000001;
G=1/50;
a1_qnl=.1e-3;a3_qnl=0; %qnl and inl are modeled using 3rd order taylor series
a1_inl=2;a3_inl=0;
Io=5; 
fs=80e3; f=1000; w1=2*pi*f; ts=1/fs; T=1/f;
h=1/(fs); K=T/ts; 

%Functions definition
syms v0 v0_mem1 t_source f
qnl_v0=a1_qnl*v0+a3_qnl*v0.^3;
qnl_v0_mem=a1_qnl*v0_mem1+a3_qnl*v0_mem1.^3;
is=Io*sin(w1*t_source);
inl_v0=a1_inl*v0+a3_inl*v0.^3;
f0=(h*G*v0)+qnl_v0+(h*inl_v0)-(h*is)-qnl_v0_mem;

% Creates the System of equations for the different tk's
for tk=1:K-1
    f(tk)=subs(f0,t_source,tk*ts);    
end

tk=K;
f(tk)=subs(f0,t_source,0);

%Initial solutions
n=1;
sol(1,:)=ones(1,K);

error=1;
jacobian=zeros(K,K);
while (abs(error)>raphson_error)
    % Constructs the jacobian
    i=1;
    for tk=1:K-1                   
        jacobian(tk,i)=subs(diff(f(tk),v0_mem1),{v0_mem1,v0},{sol(n,tk),sol(n,tk+1)});
        jacobian(tk,i+1)=subs(diff(f(tk),v0),{v0_mem1,v0},{sol(n,tk),sol(n,tk+1)});
        i=i+1;
    end
    tk=K;i=1;
    jacobian(tk,i)=subs(diff(f(tk),v0),{v0_mem1,v0},{sol(n,tk),sol(n,1)});
    jacobian(tk,tk)=subs(diff(f(tk),v0_mem1),{v0_mem1,v0},{sol(n,tk),sol(n,1)});
    inv_jacobian=inv(jacobian);
    
    % Finds F(v0)^i
    j=1;
    for tk=1:K
        if (tk==K) 
            fv0(tk)=subs(f(tk),{v0_mem1,v0},{sol(n,K),sol(n,1)});
        else 
            fv0(tk)=subs(f(tk),{v0_mem1,v0},{sol(n,j),sol(n,j+1)});
            j=j+1;
        end
    end
    
    % Intermediate solution for newton raphson
    sol(n+1,:)=sol(n,:)'-(inv_jacobian*fv0')
    error=sum(abs(sol(n+1,:)-sol(n,:)))
    n=n+1;
end

% The solution
n=n-1;
vout=[sol(n,:) sol(n,1)]; % The element K is equal to the first element

% The Source
n=1; tsource=0:ts:T;
for tk=0:ts:T
    source(n)=subs(is,t_source,tk);
    n=n+1;
end

plot(tsource,source,'-*g',tsource,vout,'-*b');
legend('source','vout');
grid on;