ODE coupled with classic equation

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

[darova]

Do you have source formula/equation?

[MartinM]

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

Answer 1

[darova]

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)

[MartinM]

1: question : The value 1 at the end should be rho0?

2 : in the first code i right I creat E linked to the time vector t and add your solution

I doubled Letter E==> EE etc

clc
clear all
close all
c=3e8;
T0=100e-15;
lambda0=515e-9;
w0=2.*pi.*c./lambda0;
Ej=100e-6;
Pp=Ej./T0;
dt=T0./100;
t=-T0*5:dt:T0*5;

w=2.*pi.*c./lambda;

E=Pp.*exp(-(t./T0).^2).*cos(w0.*t);
% plot(t,E)

a=9.9992e28;
b=3.5482e11;
rho0=2.7e26;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
EE = @(tt) Pp*exp(-tt^2/T0^2)*cos(w0*tt);
ww = @(tt) a/EE(tt)*exp(-b/EE(tt));
rho = @(tt,rho) ww(tt)*(rho-rho0);
[tt,r] = ode45(rho,[-T0*5:dt:T0*5],rho0);
plot(tt,r)

fplot(EE,tt)

Problem 1 : r(1)= rho0, but the other r = NaN

Problem 2 : EE is not found when I plot is limited xaxis at [-500e-15 -499e-15] : Solved, the interval contains to much point...don't understand the problem.

[darova]

1: question : The value 1 at the end should be rho0?

• Value at the end is initial condition for rho (not rho0) and it can't be rho=rho0, because rho-rho0=0

Problem 1 : r(1)= rho0, but the other r = NaN

Problem 2 : EE is not found when I plot is limited xaxis at [-500e-15 -499e-15]

• Try to analyze your function because there are many values where it doesn't exist (NaN of Inf)

• Here is rho function

clc,clear
c = 3e8;
T0 = 100e-15;
lambda0 = 515e-9;
w0=2.*pi.*c./lambda0;
Ej=100e-6;
Pp=Ej./T0;
dt=T0./100;
a=9.9992e28;
b=3.5482e11;
rho0=2.7e26;
EE = @(tt) Pp*exp(-tt.^2/T0^2).*cos(w0*tt);
ww = @(tt) a./EE(tt).*exp(-b./EE(tt));
rho = @(tt,rho) ww(tt)*(rho-rho0);
% [tt,r] = ode45(rho,[-T0 T0*5],rho0*1.1);
tt = linspace(-5*T0,5*T0);
rho(tt,1)'
plot(tt,rho(tt,1))

[MartinM]

Sory again,

on the last programme you provde, rho is not drho/dt ? how can i solve it?

[darova]

rho IS drho/dt

But solution cannot be found if its derivative (drho/dt) NaN or Inf

[MartinM]

EE = @(tt) Pp*exp(-tt.^2/T0^2).*cos(w0*tt);
ww = @(tt) a./EE(tt).*exp(-b./EE(tt));
rho = @(tt,rho) ww(tt)*(rho-rho0);
% [tt,r] = ode45(rho,[-T0 T0*5],rho0*1.1);
tt = linspace(-5*T0,5*T0);
rho(tt,1)'
plot(tt,rho(tt,1))

It's not r wich is drho/dt (the solution if their is)?

[darova]

• 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]

ok it's clear now,

THANKS :)

Answer 2

[MartinM]

Hi again

I update the problem. Nowtheir is an other differential euation linked to the other...How can i do this?

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

[darova]

Don't understand what are you asking

[MartinM]

How can I add this differential equation to the other ? Because it’s linked to rho and E again

I need to use the same method?

[darova]

Yes, the same method

[MartinM]

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')

[darova]