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

4 次查看（过去 30 天）
Rounak Narde
Rounak Narde 2018-8-13
编辑： Rounak Narde 2018-8-13
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!

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2018-8-13
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.
  1 个评论
Rounak Narde
Rounak Narde 2018-8-13
编辑：Rounak Narde 2018-8-13
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?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2017a

Community Treasure Hunt

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

Start Hunting!

Translated by