How to solve this differential equation inside derivation are there

1 view (last 30 days)
function kk1
a0=10; r0=70;
%main function to solve
[t,r]=ode45(@fn,[0 1200],[0.01 0.01 0.01 0.01 0.01 0.01])
plot(r(:,1),r(:,3))
function dr = fn(t,r)
% relation with actual variables
x=r(1);y=r(2); z=r(3);px=r(4);py=r(5); pz=r(6);
ta = 200;
zr=r0^2/2;
rho=sqrt((x^2+y^2)/r0^2);
eps=r0/zr;
R=z* (1+ (r0^2/(2* z))^2);
f=sqrt(1+(2*z/r0^2)^2);
l=1;
fr=(a0/f)*(sqrt(2)*rho/f)^l*exp(-rho^2/f^2)*exp(-((t-z)/ta)^2);
ph=(t-z)+2*atan(2*(z)/r0^2)-(x^2+y^2)/(2* R));
ax1=fr*cos(ph);
ay1=0;
az1=-d/dx(ax1);
bx1=0;
by1=ax1;
bz1=- d/dy(ax1);
g=sqrt(1+px^2+py^2+pz^2);
dr=zeros(6,1);
dr(1)=px/g; dr(2)=py/g; dr(3)=pz/g;
dr(4)=-ax1*(1 - pz/g)-py/g* bz1;
dr(5)=-ay1+pz/g* ay1+px/g* bz1;
dr(6)=-az1-px/g* ax1-py/g *ay1;
end
end
here bz1=is the partial derivation w r to y of ax1 similar az1 is the partial diff w r to x with in the equation so how to solve this problem

Answers (1)

Torsten
Torsten on 8 May 2022
Edited: Torsten on 8 May 2022
Execute this code separately and copy the expressions obtained for az1 and bz1 into your code from above.
If you want more variability, also declare a0, r0, l and ta as syms. You'll get an expression for az1 and bz1 depending on t,x,y,z,a0,r0,l and ta. Make az1 and bz1 functions like
function value = AZ1(t,x,y,z,a0,r0,l,ta)
value = (expression obtained from the below code)
end
function value = BZ1(t,x,y,z,a0,r0,l,ta)
value = (expression obtained from the below code)
end
and call them in your code as
az1 = AZ1(t,x,y,z,a0,r0,l,ta)
bz1 = BZ1(t,x,y,z,a0,r0,l,ta)
syms t x y z
a0=10;
r0=70;
R=z* (1+ (r0^2/(2* z))^2);
zr=r0^2/2;
eps=r0/zr;
rho=sqrt((x^2+y^2)/r0^2);
f=sqrt(1+(2*z/r0^2)^2);
ph=(t-z)+2*atan(2*(z)/r0^2)-(x^2+y^2)/(2* R);
l=1;
ta = 200;
fr=(a0/f)*(sqrt(2)*rho/f)^l*exp(-rho^2/f^2)*exp(-((t-z)/ta)^2);
ax1=fr*cos(ph);
az1 = -diff(ax1,x)
bz1 = -diff(ax1,y)
%az1 = simplify(az1)
%bz1 = simplify(bz1)
az1 = matlabFunction(az1)
bz1 = matlabFunction(bz1)

Tags

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!