[SciPy-User] port a PDE system from matab to scipy

[Massimo Di Stefano]

Hello All,


i'm tring to port a matlab code to run under scipy.


the code solve a pde system using ode45,

the matlab code is like : 

tmax = 100 ;
x_d = 60 ; x_s = 0 ;
n = 10 ;
C_0 = 10 ; V_0 = 5 ;
a = 0.9 ; p = 0.15 ; b = 0.65 ; r = 0.1 ;
x=linspace(x_s,x_d,n+2); 
h=(x_d-x_s)/(n+1);
C_x_s = 0 ; C_x_d = 0 ; V_x_s = 0 ; V_x_d = 0 ;

CViniz=zeros(2*n,1); 
CViniz(1)=C_0; 
CViniz(n+1)=V_0;

[t,CV] = ode45('AvvDiffPredRatevar',[0 tmax],CViniz); 


% where AvvDiffPredRatevar is :

function CVprimo = AvvDiffPredRatevar(t,CV)
global n h x C_x_s C_x_d V_x_s V_x_d a b p r
CVprimo(1:n) = -0.2/h*(CV(1:n)-[C_x_s; CV(1:(n-1))])+0.2/(h^2)*....
         ([C_x_s; CV(1:n-1)]-2*CV(1:n)+[CV(2:n); C_x_d])+a*CV(1:n)-b*CV(1:n).*CV((n+1):end);
CVprimo(n+1:2*n) = -0.2/h*(CV(n+1:end)-[V_x_s; CV(n+1:end-1)])+0.2/(h^2)*....
         ([V_x_s; CV(n+1:end-1)]-2*CV(n+1:end)+[CV(n+2:end); V_x_d])+p*CV(n+1:end)-r*CV(1:n).*CV(n+1:end);
CVprimo=CVprimo';

in scipy i start to write :

tmax = 100 ;
x_d = 60 ;
x_s = 0 ;
n = 10 ;
C_0 = 10 ;
V_0 = 5 ;
a = 0.9 ;
p = 0.15 ;
b = 0.65 ;
r = 0.1 ;
x=linspace(x_s,x_d,n+2); 
h=(x_d-x_s)/(n+1);
C_x_s = 0 ;
C_x_d = 0 ;
V_x_s = 0 ;
V_x_d = 0 ;

CViniz=zeros(2*n,1); 
CViniz(1)=C_0; 
CViniz(n+1)=V_0;


scip now i have to define the "ratefunction" AvvDiffPredRatevar(t,CV) 
have you any suggestion on how to translate the matlab ratefunciotn code to python ?
i suppose i have to write a new python function : def AvvDiffPredRatevar(t,CV)
and then i have to pass it to scipy : 
scipy.integrate.ode(f).set_integrator('dopri5') where "f" is "AvvDiffPredRatevar(t,CV)"
make sense for you ? or i'm wrong ?

thanks for any hints!!!

Best regards,

Massimo.





-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100928/7b9e9cde/attachment.html>