Inverse Problem
This demo is implemented in a single Python file. Download here: tutorial_inverse.ipynb
This demo illustrates how to:
How to use
InverseProblem
.:from module1_mesh import* from module2_forward import* from module3_inverse import* from module4_auxiliar import* import matplotlib.pyplot as plt
Mesh
mesh_inverse, mesh_direct=MyMesh(r=1, n=8, n_vertex=321)
mesh_direct=refine(mesh_direct)
plt.figure(figsize=(8, 8))
plt.subplot(1,2,1)
plot(mesh_direct);
plt.subplot(1,2,2)
plot(mesh_inverse);
![../../_images/mesh10.png](../../_images/mesh10.png)
Current
"Current"
n_g=6
I_all=current_method(n_g, value=1, method=1)
"Checking hyphotesis"
Verifyg(I_all, mesh_direct)
Verifyg(I_all, mesh_inverse)
Defining Gamma function
"Defining gamma"
ValuesCells0=GammaCircle(mesh_direct,3,1,0.35, 0.35, 0.30); #malha, cond_interna, cond_externa, raio.
"Plot"
V_DG=FiniteElement('DG',mesh_inverse.ufl_cell(),0)
gamma_direct=plot_figure(mesh_direct, V_DG, gamma0, name="Resposta gamma");
![../../_images/gamma10.png](../../_images/gamma10.png)
Forward Problem
"ForwardProblem"
VD=FiniteElement('CG',mesh_direct.ufl_cell(),1)
ForwardObject=ForwardProblem(mesh_direct)
list_u0=ForwardObject.solve_forward(VD, gamma0, I_all)
u0_boundary=ForwardObject.boundary_array(mesh_inverse) #Boundary data in the mesh_inverse
Noise plot
"Noise Parameters"
noise_level=0.005 #0.5%
noise_type='uniform'
seed=1
u0_boundary=ForwardObject.add_noise(noise_level, noise_type, seed, mesh_inverse)
end=len(u0_boundary)
if noise_level>1E-5:
u0_boundary_c=ForwardObject.boundary_array(mesh_inverse, concatenate=False)
u_dados_noise=u0_boundary[0:int(end/n_g)]
data_u0=plot_boundary(mesh_inverse, data=u0_boundary_c[index], name='boundary u0_0', line=0, plot=False)
data_u0_noise=plot_boundary(mesh_inverse, data=u_dados_noise, name='boundary u0_0', line=0, plot=False)
plt.figure(figsize=(14, 5))
plt.subplot(1,2,1)
plt.plot(data_u0_noise[:,0], data_u0_noise[:,1], linewidth=0, marker='.', markersize=2, label="noise", color='red');
plt.plot(data_u0[:,0], data_u0[:,1], linewidth=0, marker='.', markersize=2, label="no_noise", color="blue");
plt.legend()
plt.subplot(1,2,2)
plt.plot(data_u0[:,0], data_u0[:,1]-data_u0_noise[:,1], linewidth=1, marker='.', markersize=2, label="noise", color="blue");
plt.legend()
plt.show()
![../../_images/noise3.png](../../_images/noise3.png)
Inverse Problem
Inverse Object and Parameters
"Object"
VI=FiniteElement('CG',mesh_inverse.ufl_cell(),1)
InverseObject=InverseProblem(mesh_inverse, VI, u0_boundary, I_all)
"Solver Parameters"
InverseObject.set_solverconfig(step_limit=100)
InverseObject.set_NewtonParameters(mu_i=0.97, mumax=0.999, nu=0.99, R=0.999)
#InverseObject.set_InnerParameters(inner_method='Landweber', land_a=20)
InverseObject.set_InnerParameters(inner_method='ME', ME_reg=5E-3)
#InverseObject.set_InnerParameters(inner_method='CG')
#InverseObject.set_InnerParameters(inner_method='LM', LM_c0=0.95, LM_q=0.05)
#InverseObject.set_InnerParameters(inner_method='Tikhonov', Tik_c0=0.95, Tik_q=0.05)
"Noise Parameters"
tau=1.2
InverseObject.set_NoiseParameters(tau=tau, noise_level=0.005)
print(f"Discrepance stop {tau*noise_level*100}% residue")
"Extra Parameters"
InverseObject.set_answer(gamma0, mesh_direct)
InverseObject.weight_value=False
#MyCellsgamma_k=np.ones(mesh_inverse.num_cells())*2.0
#InverseObject.set_firstguess(MyCellsgamma_k)
Discrepance stop 0.6% residue
#Solver
InverseObject.solve_inverse()
Getting Results
"Getting the result"
gamma_k=InverseObject.gamma_k
"Getting data_summary"
innerstep_vec=InverseObject.innerstep_vec
res_vec=InverseObject.res_vec
mun_vec=InverseObject.mun_vec
error_vec=InverseObject.error_vec
"Plot result"
V_DG=FiniteElement('DG',mesh_inverse.ufl_cell(),0)
gamma_k=plot_figure(mesh_inverse, V_DG, gamma_k, name=" ");
![../../_images/gammares1.png](../../_images/gammares1.png)
Fancy plot
"Fancy plot"
cmap='jet'
#cmap='inferno'
fig=plt.figure(figsize=(9, 9))
plt.subplot(1,2,1);
p=plot(gamma_direct);
p.set_cmap(cmap)
plt.colorbar(p,fraction=0.046, pad=0.0);
plt.subplot(1,2,2);
p=plot(gamma_k);
p.set_cmap(cmap);
plt.colorbar(p,fraction=0.046, pad=0.0);
fig.tight_layout()
![../../_images/fancy.png](../../_images/fancy.png)
Plot error and residue
plt.figure(figsize=(12, 4));
"Subplot 1"
plt.subplot(1,2,1)
plt.ylabel('Residue (%)');
plt.xlabel('Step');
plt.title('Residue');
plt.plot(range(0,len(res_vec)), res_vec);
"Subplot 2"
plt.subplot(1,2,2)
plt.ylabel('Error(gamma_k - gamma_exata) (%)');
plt.xlabel('Step');
plt.title('Error');
plt.plot(range(0,len(error_vec)), error_vec);
![../../_images/reserror.png](../../_images/reserror.png)
Summary data
summary_data=np.vstack((error_vec, res_vec, innerstep_vec, mun_vec))
data=pd.DataFrame(summary_data.T)
experiments=["Error (%)", "Residue (%) ", "Inner Step", "mu_n"]
data.columns = experiments
data
![../../_images/summary1.png](../../_images/summary1.png)