[pydstool-users] help with nonlinear DAE

[<fre...@gm...>]

Hya,

currently investigateing PyDSTool to solve some DAEs. I tried to implement an easy electric circuit with a capacitor discharging through a non-linear resistor. (I know that this system could be reduced to one ODEs but I'm practising.) Below the code for PyDSTool and also for matlab. Both work for a linear resistor, but PyDSTools ceases to work for even tiny non-linearities whereas matlab works fine.

where am I going wrong? Any options I need to set? 

Thanks for the help

mauro
--------------------------------------------

#! /usr/local/bin/python
"""Simple example of solving a differential-algebraic equation.
A capacitor and non-linear  resistor in series.

Resistor:
- delta_v = R*I^non_lin*sign(I)
Capacitor:
I = C * d(delta_v)/dt

for non_lin=1 has solution:
delta_v = exp(t/(-RC)) * delta_v_0
I = -1/R( exp(t/(-RC)) * delta_v_0)
"""

from PyDSTool import *

DSargs = args()

# parameters
R = 80.                                 # resistance
C = 1e-3                                # capcitance
non_lin = 1.01                          # non linearity
DSargs['pars'] = {'R': R,
                  'C': C,
                  'non_lin': non_lin}

# IC: making sure they satisfy the algebraic constraint
delta_v_0 = 100
I_0 = -sign(delta_v_0)*(delta_v_0/R)**(1./non_lin)
DSargs['ics'] = {'I': I_0, 'delta_v': delta_v_0}


# eqn. definition
DSargs['varspecs'] = {'I':
                      'I*R*pow(abs(I),non_lin-1) + delta_v',
                      'delta_v':
                      'I/C'}
DSargs['fnspecs'] = {'massMatrix':
                     (['t','I','delta_v'], '[[0,0],[0,1]]')}
DSargs['vars'] = ['I', 'delta_v']

# numerical parameters
DSargs['algparams'] = {'init_step': 0.0005, 'max_step': 0.1,
                       'rtol': 1e-9, 'atol': 1e-9,
                       'max_pts': 1000000}
DSargs['checklevel'] = 1
DSargs['name'] = 'DAE_nonlin_resistance_capacitor'

dae = Generator.Radau_ODEsystem(DSargs)
dae.set(tdata=[0, 1])
traj = dae.compute('nonlin')
pd = traj.sample(dt=0.05)

# exact solutions for non_lin=1
I = lambda t : -1/R*( exp(t/(-R*C)) * delta_v_0)
delta_v = lambda t : exp(t/(-R*C)) * delta_v_0

figure()
plot(pd['t'], pd['I'])
hold(True)
plot(pd['t'], I(pd['t']), 'bd')
twinx()
plot(pd['t'], pd['delta_v'], 'g')
plot(pd['t'], delta_v(pd['t']), 'gd')
show()

----------------------------------------------------

function [t,out] = resitor_capacitor();
% to integrate the system of DAE arrising from a non 
% linear resistor in series with a capacitor

plot_yes = 1;  

R = 80;
C = 1e-3;
non_lin = 1.01;

%IC
delta_v_0 = 100
I_0 = -sign(delta_v_0)*(delta_v_0/R)^(1/non_lin)

ic = [I_0 delta_v_0]';

tspan = [0 1]';    % integration time
M = diag([0 1]);    % mass matrix

opts = odeset('Mass',M,'MassSingular','yes', ...
              'stats','on','MStateDependence','none');
[t,out] = ode15s(@int_fun,tspan,ic,opts);

if plot_yes == 1
    figure
    subplot(211)
    plot(t, out(:,1))
    subplot(212)
    plot(t, out(:,2),'g')
end
    
% function to be integrated:
    function out = int_fun(t,in)
        % assign variables for convinience
        I = in(1);
        delta_v = in(2);


        % do the maths:
        out(1) = I*R*(abs(I))^(non_lin-1) + delta_v;
        out(2) = I/C;
        
        out = out';
       
    end

end
-- 
GMX startet ShortView.de. Hier findest Du Leute mit Deinen Interessen!
Jetzt dabei sein: http://www.shortview.de/?mc=sv_ext_mf@gmx