ODE coupled with classic equation

1 次查看（过去 30 天）
MartinM
MartinM 2019-10-23
评论： darova 2019-10-23
Hi everybody.
After some research can't found a solution..
I have 2 variable wich depend on time : E and W(E)
then I have an differential equation of rho inked to W so linked to E so linked to t.
Can I use E et W as vector inside the ODE declaration?
clc
clear all
close all
u=2.405;
c=3e8;
T0=100e-15;
lambda0=515e-9;
w0=2.*pi.*c./lambda0;
Ej=100e-6;
Pp=Ej./T0;
r=18e-6;
th=250e-9;
s=0.085;
dt=T0./1000;
t=-T0*5:dt:T0*5;
Fs=1./dt;
nn=length(t),
freq = Fs*linspace(0,(nn/2),(nn/2)+2)/nn+c/lambda0;
freq=fliplr(freq(1:end-1));
l=c./freq;
ll=-fliplr(l);
lll=ll-ll(1)+l(end);
lll = (circshift(lll',-1))';
lambda=[l lll];
lambda=lambda(1:end-1);
w=2.*pi.*c./lambda;
E=Pp.*exp(-(t./T0).^2).*cos(w0.*t);
% plot(t,E)
a=r.*( 1+ s.*(2*pi.*c).^2./ (w.*w.*r.*th) ).^(-1);
% plot(lambda,a)
a=9.9992e28;
b=3.5482e11;
rho0=2.7e26;
W=a./(abs(E)).*exp(-b./(abs(E)));
syms rho(t) EE(t) WW(t)% Y ;
ode1= EE== Pp.*exp(-(t./T0).^2).*cos(w0.*t);
ode2 = WW==a./(abs(EE)).*exp(-b./(abs(EE)))
ode3 = diff(rho,t) == W(t) .*(rho0 - rho);
ode=[ode1 ode2] ode3
rhoSol=solve(ode)
%%%or
yms rho(t) ;
ode = diff(rho,t) == W(t) .*(rho0-rho);
rhoSol=solve(ode)
If you have an idea to solve this?
Regards
MM
  2 个评论
darova
darova 2019-10-23
Do you have source formula/equation?
MartinM
MartinM 2019-10-23
编辑：MartinM 2019-10-23
sure , t is my time vector from the code above, tau et w0 are initial values from the code aboce
Hope that can help.
I think it's simple, but W is a vector from the vector E.
Regards

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

darova
darova 2019-10-23
Try this
E = @(t) E0*exp(-t^2/tau^2)*cos(w0*t);
w = @(t) a/E(t)*exp(-b/E(t));
rho = @(t,rho) w(t)*(rho-rho0);
[t,r] = ode45(rho,[0 0.1],1);
plot(t,r)
  7 个评论
显示 5更早的评论
darova
darova 2019-10-23
  • It's not r wich is drho/dt (the solution if their is)?
Yes. BUt if rho is inf or NaN drho/dt cannot be found. You asked why all r are NaN - this is the answer, because of rho
MartinM
MartinM 2019-10-23
ok it's clear now,
THANKS :)

请先登录，再进行评论。

更多回答（1 个）

MartinM
MartinM 2019-10-23
Hi again
I update the problem. Nowtheir is an other differential euation linked to the other...How can i do this?
tout.png
clc
clear all
close all
u=2.405;
c=3e8;
%%%LASER%%%
T0=5e-15;
lambda0=515e-9;
w0=2.*pi.*c./lambda0;
Ej=100e-6;
Pp=Ej./T0;
%%%Fibre%%%
r=18e-6;
th=250e-9;
s=0.085;
%%%GAZ%%%
Pg=10;
B1=2.6102828e-4
B2=5.694682e-4;
C1=2.01e-6;
C2=1.0043e-2;
Khi3=2.2e-25;
Epsi0=8.85418782e-12;
%%%Espace%%
dt=T0./1000;
t=-T0*5:dt:T0*5;
Fs=1./dt;
nn=length(t),
freq = Fs*linspace(0,(nn/2),(nn/2)+2)/nn+c/lambda0;
freq=fliplr(freq(1:end-1));
l=c./freq;
ll=-fliplr(l);
lll=ll-ll(1)+l(end);
lll = (circshift(lll',-1))';
lambda=[l lll];
lambda=lambda(1:end-1);
w=2.*pi.*c./lambda;
%%%Champ
E=Pp.*exp(-(t./T0).^2).*cos(w0.*t);
% plot(t,E)
%%%RAYON%%%
a=r.*( 1+ s.*(2*pi.*c).^2./ (w.*w.*r.*th) ).^(-1);
% plot(lambda,a)
lmicro=lambda.*1e6;
ng=sqrt( 1+ Pg.*( (B1.* lmicro.^2./(lmicro.^2 -C1))+ (B2.* lmicro.^2./(lmicro.^2 -C2)) ) );
betaw=( (w.*ng./c).^2 - (u./a).^2);
P=Epsi0.*Khi3.*(abs(E)).^2 .*E;
alpha=9.9992e28;
b=3.5482e11;
rho0=2.7e26;
W=alpha./(abs(E)).*exp(-b./(abs(E)));
EE = @(tt) Pp*exp(-tt.^2/T0^2).*cos(w0*tt);
ww = @(tt) alpha./abs(EE(tt)).*exp(-b./abs(EE(tt)));
rho = @(tt,rho) ww(tt)*(rho-rho0);
[tt,RHO] = ode45(rho,t,0);
plot(t,RHO)
Itry to use syms, woth somethong like...
syms RRHO(tt) EE(tt) WW(tt) JJ(tt)
ode1 = EE== Pp*exp(-tt.^2/T0^2).*cos(w0*tt);
ode2 = WW== alpha./abs(EE).*exp(-b./abs(EE));
ode3 = diff(YY) == (YY);
ode=......
%%%%not working
problem with differential and classical I think
  5 个评论
显示 3更早的评论
MartinM
MartinM 2019-10-23
Again not working.
I think the problem comes from how jj is written.Matlab would like an analyticexpression instead of a the vector RHO. I tried to change RHO by rho(tt)...
clc
clear all
close all
u=2.405;
c=3e8;
q=1.60217662e-19;
m=9.109e-31;
%%%LASER%%%
T0=5e-15;
lambda0=515e-9;
w0=2.*pi.*c./lambda0;
Ej=100e-6;
Pp=Ej./T0;
%%%Fibre%%%
r=18e-6;
th=250e-9;
s=0.085;
%%%GAZ%%%
Pg=10;
B1=2.6102828e-4
B2=5.694682e-4;
C1=2.01e-6;
C2=1.0043e-2;
Khi3=2.2e-25;
Epsi0=8.85418782e-12;
%%%Espace%%
dt=T0./1000;
t=-T0*5:dt:T0*5;
Fs=1./dt;
nn=length(t),
freq = Fs*linspace(0,(nn/2),(nn/2)+2)/nn+c/lambda0;
freq=fliplr(freq(1:end-1));
l=c./freq;
ll=-fliplr(l);
lll=ll-ll(1)+l(end);
lll = (circshift(lll',-1))';
lambda=[l lll];
lambda=lambda(1:end-1);
w=2.*pi.*c./lambda;
%%%Champ
E=Pp.*exp(-(t./T0).^2).*cos(w0.*t);
% plot(t,E)
%%%RAYON%%%
a=r.*( 1+ s.*(2*pi.*c).^2./ (w.*w.*r.*th) ).^(-1);
% plot(lambda,a)
lmicro=lambda.*1e6;
ng=sqrt( 1+ Pg.*( (B1.* lmicro.^2./(lmicro.^2 -C1))+ (B2.* lmicro.^2./(lmicro.^2 -C2)) ) );
betaw=( (w.*ng./c).^2 - (u./a).^2);
P=Epsi0.*Khi3.*(abs(E)).^2 .*E;
alpha=9.9992e28;
b=3.5482e11;
rho0=2.7e26;
W=alpha./(abs(E)).*exp(-b./(abs(E)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
EE = @(tt) Pp*exp(-tt.^2/T0^2).*cos(w0*tt);
ww = @(tt) alpha./abs(EE(tt)).*exp(-b./abs(EE(tt)));
rho = @(tt,rho) ww(tt)*(rho0-rho);
[tt,RHO] = ode45(rho,t,0);
RHO=RHO;
return
jj = @(tt,jj) q.*q./m.*RHO.*EE(tt) ;%-jj./Tc;
[tt,J] = ode45(jj,t,0);
return
plot(t,RHO,'color','r')

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by