Asked by Rounak Narde
on 13 Aug 2018

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!

Answer by 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.

Rounak Narde
on 13 Aug 2018

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?

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## KSSV (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/414609-solving-2nd-order-linear-ode-wave-equation-the-solution-at-final-values-overshoots-diverging-fro#comment_599048

## Rounak Narde (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/414609-solving-2nd-order-linear-ode-wave-equation-the-solution-at-final-values-overshoots-diverging-fro#comment_599224

Sign in to comment.