ODE coupled with classic equation

2019 10 23
2 个回答
回答已采纳
更新时间：2019 10 23
13 次查看（30 天）

0 个投票

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

0 个投票

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更早的评论 隐藏 5更早的评论

MartinM
MartinM 2019-10-23
编辑：MartinM 2019-10-23
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
darova 2019-10-23
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)
  • img1.png 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
MartinM 2019-10-23
Sory again,
on the last programme you provde, rho is not drho/dt ? how can i solve it?
darova
darova 2019-10-23
rho IS drho/dt
But solution cannot be found if its derivative (drho/dt) NaN or Inf
MartinM
MartinM 2019-10-23
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
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

0 个投票

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更早的评论 隐藏 3更早的评论

darova
darova 2019-10-23
Don't understand what are you asking
MartinM
MartinM 2019-10-23
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
darova 2019-10-23
Yes, the same method
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')

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by