Solving 2nd order linear ode - wave equation - the solution at final values overshoots diverging from theoretical solution

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!

[KSSV]

Can you show us the governing equation used for ODE45.

[Rounak Narde]

where freq, mu_0 and sigma_bulk is constant. B is magnetic field and x is depth inside the conductor.

Answer 1

[Torsten]

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.

[Rounak Narde]

Thanks Torsten for your reply!

The error does not go away. If you will increase the range of the depth (dmax) more than 20um (let's say 80um), you will start seeing the overshoot again.

This code is just a part of my actual work. My work requires multiple frequency solutions upto 100GHz. As freq go near to 100GHz or even 20GHz, the overshoots overshadows the required solution, and shoots above 10^40.

I have just a guess. can this overshoot be because of very very low value of magnetic field, when it starts to overshoot the theoretical value goes below 10^-16 or so? ode45 is unable to go that low?