Code not giving correct result when constants are given of the order of 10^-10, whilst when i give constant of the order of 10, it runs.

조회 수: 1 (최근 30일)
Following is the code for solving coupled diff eqn using rk4, it is for pupulation dynamics of a two level laser system. Here the real values of h=0.02e-9, xn=2e-9, A=5e-15, B=2e-9. With these the code does not furnish proper result, but with the following valued, it runs.
Also i don't know how to use vector function handle hence the lengthy version.
clc
clf
syms h xn A B ;
Inputh=(0.02);
Inputxn=(2);
InputA=(5);
InputB=(2);
h=vpa(Inputh)
h = 
0.02
xi=0;
xn=vpa(Inputxn);
x=vpa(linspace(xi,xn,100));
A=vpa(InputA);
B=vpa(InputB);
f_1=@(y,z)(-A*y+A*z+z/B)
f_1 = function_handle with value:
@(y,z)(-A*y+A*z+z/B)
f_2=@(y,z)(A*y-A*z-z/B)
f_2 = function_handle with value:
@(y,z)(A*y-A*z-z/B)
y=zeros(1,length(x));
z=zeros(1,length(x));
y(1)=100;
z(1)=0;
for i=1:(length(x)-1)
k_1=h*(f_1(y(i),z(i)));
l_1=h*(f_2(y(i),z(i)));
k_2=h*(f_1(y(i)+k_1/2,z(i)+l_1/2));
l_2=h*(f_2(y(i)+k_1/2,z(i)+l_1/2));
k_3=h*(f_1(y(i)+k_2/2,z(i)+l_2/2));
l_3=h*(f_2(y(i)+k_2/2,z(i)+l_2/2));
k_4=h*(f_1(y(i)+k_3,z(i)+l_3));
l_4=h*(f_2(y(i)+k_3,z(i)+l_3));
K=(k_1+2.*k_2+2.*k_3+k_4)./6;
L=(l_1+2.*l_2+2.*l_3+l_4)./6;
y(i+1)=y(i)+K;
z(i+1)=z(i)+L;
end
hold on
plot(x,y)
plot(x,z)
hold off

답변 (1개)

Sandeep
Sandeep 2023년 3월 2일
Hi Vipul kumar,
It is my understanding that you could obtain the Graph given by using the values without exponential notation but the same graph is not produced when the values used are in the order of 10^(-9).
This discrepancy is due to the fact that both the graphs use the same X-limit and Y-limit and points plotted in the order of 10^(-9) become insignificant when the Y-axis limit is in the order of 10. You can refer the following documentation page to get more insight on setting custom X and Y limits.

카테고리

Help CenterFile 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