Solving 2nd order linear ode - wave equation - the solution at final values overshoots diverging from theoretical solution
2 views (last 30 days)
Show older comments
Hi I am trying to solve second order ode - wave equation - using ode45.
where B is magnetic field; x is depth inside the conductor freq, mu_0 and sigma_bulk are constants.
Solved (theoretically):
Below is the code
% % Solving for magnetic field at the interface of conductor and air
% %
close all;
clear all;
%
freq=[15]*1E9; % Frequency 15GHz
d=[0:0.1:20]*1E-6; % depth
%
% Constants
mu_0=4*pi*1E-7; % Free space permeability
sigma_bulk=5.8E7; % conductivity of copper
%
f1=figure;
ax1=axes;
%
del=sqrt(2/(2*pi*freq*mu_0*sigma_bulk)); % Skin depth
%initial condition - Normalized Magnetic field so initial_condition(1)=1
% initial_condition(2) is found theoretically from theoretical solution.
initial_condition=[1;-(1+1i)/del];
%
% ODE - to be solved numerically
diff_func_smooth = @(x,B) [B(2);1i*2*pi*freq*mu_0*sigma_bulk*B(1)];
%
% Smooth Surfaces
% options = odeset('AbsTol',[1e-16 1e-16],'RelTol',1e-10, 'MaxStep',1E-5);
% [t_smooth,B_smooth]=ode45(diff_func_smooth, [0 20]*1E-6, initial_condition,options);
[t_smooth,B_smooth]=ode45(diff_func_smooth, d, initial_condition);
plot(ax1,t_smooth,abs(B_smooth(:,1)),'-k')
% plot(ax1,t_smooth,abs(B_smooth(:,1)),'or')
hold(ax1,'on')
grid(ax1,'on');
%
% Theoretical Normalized Magnetic field
b_smooth_solved=exp(-(1+1i)*d/del);
plot(ax1,d,abs(b_smooth_solved),'--g')
ylabel(ax1, 'Normalized B')
xlabel(ax1, 'Depth')
legend(ax1, 'Solved by ode45', 'Theoretical');
The final plot is
At the final depth values around 20um. The magnetic field overshoots despite the theoretical normalized magnetic field attending towards 0.
I was trying to use options, RelTol, AbsTol, MaxStep. Also, I have tried different solvers ode23, ode23s. The error at the end increases and decreases. But error remains. I am not sure which option to use, or what I am doing wrong. Please help Thanks!
2 Comments
Answers (1)
Torsten
on 13 Aug 2018
Numerical problem.
If you choose the product "2*pi*freq*mu_0*sigma_bulk" to be one power of ten less (e.g freq = 15e8 instead of freq = 15e9), all works well.
Best wishes
Torsten.
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!