## pydstool-users

 [pydstool-users] help with nonlinear DAE From: - 2008-05-20 14:21:05 ```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@... ```