Solving Second Order Differential Equation--Getting Nonsensical Answers!

1 次查看（过去 30 天）
J Shaw
J Shaw 2012-5-24
Hi Everyone!
I am trying to solve the following equation using Matlab:
d^2(r)/dz^2 + (1/gamma)*dr/dz*dgamma/dz + (constant/gamma)*r = 0.
where gamma is a function of z but dgamma/dz is constant. I have set it up as shown below. However, my result appears to not depend on the dr/dz term. It gives the same result with or without that term! Any ideas as to what I am doing wrong?
Thanks so much!
~~~~~Code to establish function~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%===WRITE EQUATION AS FIRST ORDER SYSTEM===================================
%Let r1(z) be r(z)
%Let r2(z) be dr(z)/dz
%Let r3(z) be gamma
%Therefore r1prime(z)=dr(z)/dz=r2(z)
%Therefore our system of equations becomes
%==>r1prime = r' = r2
%==>r2prime = r'' = -(gamma^-1)*r2*dgamma-((w_p^2)/(2*gamma*c^2))*r1
%==>r3prime = r_3' = gamma-gamma_0
%===DEFINE AN M-FILE WITH THE SYSTEM OF EQUATIONS==========================
function rprime=M12_05_24_kBetaFunction(z,r)
%===USER INPUT=============================================================
plasma_density=3e19;
z_first_limit=500e-
lambda_0=800e-9;
%Choose ramp profile
n_ramp=0; %Density step
%===DEFINE CONSTANTS=======================================================
e=1.6e-19; %Charge of electron [C].
c=3e8; %Speed of light [m/s].
epsilon=8.85e-12; %Permittivity of free space
m=9.11e-31; %Mass of electron [kg]
%===CALCULATE VALUES=======================================================
%Convert plasma density to SI units
plasma_density=plasma_density*100^3; %Convert plasma density to [m^-3].
%Set plasma density and energy gradient for each region
if z<z_first_limit
n=plasma_density; %Set density inside plasma
EnergyGradient=1.54e-20*sqrt(n); %Set energy gradient inside plasma
%EnergyGradient=0;
else
n=n_ramp; %Sets density on ramp
EnergyGradient=0; %Assume no energy gain on ramp
end
%Define plasma frequencies
omega_0=2*pi*c/lambda_0;
omega_p=sqrt((e^2*n)/(m*epsilon));
omega_p0=sqrt((e^2*plasma_density)/(m*epsilon));
%Define gamma and energy at trapping
gamma_0=omega_0/omega_p0;
Energy_0=.511*(gamma_0-1)*1.6e-13; %[J]
%Set energy and gamma for each region
if z<z_first_limit
Energy=Energy_0+EnergyGradient*z; %Calculate energy at each z [J]
EnergyMeV=Energy/(1.6e-13); %Convert energy back to [MeV]
gamma=EnergyMeV/.511+1; %Calcualte gamma at each z
else
Energy=Energy_0+1.54e-20*sqrt(plasma_density)*z_first_limit; %Energy remains at the
%same value as when it leaves the plasma [J].
EnergyMeV=Energy/(1.6e-13); %Convert energy back to [MeV]
gamma=EnergyMeV/.511+1; %Calcualte gamma at each z
end
rprime=zeros(2,1);
rprime(1) = r(2);
rprime(2) = -(1/gamma)*r(2)*(1000/.511)*EnergyGradient-(omega_p^2)*r(1)/(2*c^2*gamma);
~~~~~Code to evaluate function~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%===CLOSE AND CLEAR EVERYTHING=============================================
close all
clear all
clc
%===USER INPUT=============================================================
lambda_0=800e-9; %Laser wavelength [m]
plasma_density=3e19; %Plasma density [cm^-3]
steps=100;
z_start=1e-6; %Location that electrons start oscillating [m]
z_exit=500e-6; %Location that plasma ends [m]
z_end=200e-6; %Length that electrons propagate in vacuum [m]
%Enter initial values of r1 and r2 (where again, r1=r(z) and r2=dr/dz).
r1_0=3e-6; %Distance off axis where electron bunch starts [m].
r2_0=0;
%===SET CONSTANTS==========================================================
e=1.6e-19; %Charge of electron [C]
c=3e8; %Speed of light [m/s].
epsilon=8.85e-12; %Permittivity of free space
m=9.11e-31; %Mass of electron [kg].
%===CALCULATE VALUES=======================================================
%Calculate plasma frequencies
omega_0=2*pi*c/lambda_0;
omega_p0=sqrt((e^2*plasma_density*100^3)/(m*epsilon));
%Define gamma at trapping
gamma_0=omega_0/omega_p0;
%Set initial values
r_0=[r1_0, r2_0]; %r_0(1) is r1_0 [m]. r_0(2) is r2_0 [dimenionsless].
ztotal=zeros(1,2*steps);
ztotal(1:steps)=linspace(z_start,z_exit,steps);
dz=(z_exit-z_end)/steps;
ztotal(steps+1:2*steps)=linspace(z_exit+dz,z_end+z_exit,steps);
%===ACTUALLY RUN ODE=======================================================
[Z,Result]=ode45(@M12_05_24_kBetaFunction,ztotal,r_0);
%Note: "Z" is all the values of z. "Result(:,1)" is all the
%values of r. "Results(:,2)" is all the values of dr/dz.

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

回答（1 个）

Walter Roberson
Walter Roberson 2012-5-24
You indicate that dgamma/dz is constant . Is that the same constant as in your equation?
If they are then I find an analytic solution of
r(z) = C2 * BesselJ(0, 2*((constant*z + C1)/constant)^(1/2)) +
C3 * BesselY(0, 2*((constant*z + C1)/constant)^(1/2))
gamma(z) = constant*z + C1
In the above, C1, C2, C3 are constants of integration, and "constant" is the constant that dgamma(z)/dz is equal to.

类别

在 Help Center 和 File Exchange 中查找有关 Plasma Physics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by