迭代法求方程的根 (ψ(x)-x=0的根)

本文介绍了一种使用迭代法求解非线性方程ψ(x)-x=0的不动点的方法,并提供了一个MATLAB函数实现。该函数接受一个函数句柄、初始猜测值和容许误差作为输入参数,返回不动点的近似值及其迭代次数。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

function [x_star,k]=my_iterate(fun,a,ep)


%用非线性求解迭代法通过 x=ψ(x),求解ψ(x)-x=0的根
%N最大迭代次数
% a是不动点的近似
%ep容许误差


N=500
my_fun=fun;
k=1;
x_0=a
x_1=feval(my_fun,a)

    while abs(x_1-x_0) && N >= 0
        k=k+1;
        N=N-1
        x_0=x_1
        x_1=feval(my_fun,x_0)
       
    end

    x_star=x_1;

end


如求解x=exp(-x)

fun=inline('exp(-x)')

[x_star,k]=my_iterate(fun,0.5)

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import fsolve from scipy.integrate import trapz from scipy.interpolate import interp1d import os # ============================================================= # 物理常数 # ============================================================= # ε0: 真空介电常数 (F/cm) epsilon_0 = 8.854e-14 # q: 电子电荷 (C) q = 1.602e-19 # k: 玻尔兹曼常数 (J/K) k = 1.380649e-23 # T: 温度 (K) T = 300 # β = q/(kT): 热电压的倒数 (V^{-1}) beta = q / (k * T) # ============================================================= # 材料参数 (据文件2.pdf更新) # ============================================================= params = { 'N_D': 1e16, # cm^{-3}, 衬底施主掺杂浓度 (公式4.4) 't_F': 200e-7, # cm, 铁电层厚度 (200 nm) 'epsilon_F': 200, # 铁电层相对介电常数 'epsilon_i': 3.9, # 绝缘层相对介电常数 (公式4.4) 'epsilon_Si': 11.7, # Si相对介电常数 #'P_s': 10e-6, # C/cm^2, 饱和极化 #'P_r': 8.5e-6, # C/cm^2, 剩余极化 #'E_c': 50e3, # V/cm, 矫顽电场 #'delta': 75e3, # V/cm, 铁电材料参数 'EOT': 4e-7, # cm, 等效氧化层厚度 'A_f': 100e-8, # cm^2, 栅面积 'W': 10e-4, # cm, 沟道宽度 'L': 1e-4, # cm, 沟道长度 'mu': 200, # cm^2/V·s, 载流子迁移率 'n_i': 1.5e10, # cm^{-3}, 本征载流子浓度 } # ============================================================= # 计算派生参数 # ============================================================= # 绝缘层电容 C_i = A_f * ε_i * ε_0 / EOT (公式4.7) params['C_i'] = params['A_f'] * epsilon_0 * params['epsilon_i'] / params['EOT'] # 德拜长度 L_D = sqrt(ε_Si ε_0 φ_t / (q N_D)) (公式4.4) phi_t = k * T / q # 热电压 params['L_D'] = np.sqrt(epsilon_0 * params['epsilon_Si'] * phi_t / (q * params['N_D'])) # 费米势 ψ_b = -φ_t * ln(N_D / n_i) (公式4.10) params['psi_b'] = -phi_t * np.log(params['N_D'] / params['n_i']) # ============================================================= # 读取外部P-E关系数据 # ============================================================= pe_file = r"D:\code\phi\pe.txt" if not os.path.exists(pe_file): raise FileNotFoundError(f"P-E关系文件未找到: {pe_file}") # 加载数据 pe_data = np.loadtxt(pe_file) E_ext = pe_data[:, 0] # MV/m P_ext = pe_data[:, 1] # C/m² # 将数据拆分为上扫和下扫分支 max_idx = np.argmax(E_ext) E_up = E_ext[:max_idx+1] # 上扫分支 P_up = P_ext[:max_idx+1] E_down = E_ext[max_idx:] # 下扫分支 P_down = P_ext[max_idx:] # 转换为代码使用的单位 E_up = E_up *28.05* 1e4 # MV/m -> V/cm P_up = P_up * 0.26* 1e-4 # C/m² -> C/cm² E_down = E_down *28.05*1e4 P_down = P_down * 0.26 *1e-4 # 创建插值函数 pe_interp_up = interp1d(E_up, P_up, kind='linear', fill_value="extrapolate") pe_interp_down = interp1d(E_down, P_down, kind='linear', fill_value="extrapolate") # ============================================================= # 外部极化函数 # ============================================================= def external_polarization(E, direction): if direction == 'up': return pe_interp_up(E) else: # down return pe_interp_down(E) # ============================================================= # C-V 特性仿真 (对应文件2.pdf第4.1节) # ============================================================= def Q_s(psi_s): """ 半导体表面电荷密度 (C/cm^2) 据文件公式(4.4): Q_s(ψ_s) = -sign(ψ_s) * [√2 * ε0 * ε_i] / (β * L_D) * sqrt[ (n_i²/N_D²)(e^{-βψ_s} + βψ_s - 1) + (e^{βψ_s} - βψ_s - 1) ] """ term = (params['n_i']**2 / params['N_D']**2) * (np.exp(-beta * psi_s) + beta * psi_s - 1) + \ (np.exp(beta * psi_s) - beta * psi_s - 1) return -np.sign(psi_s) * (np.sqrt(2) * epsilon_0 * params['epsilon_i']) / (beta * params['L_D']) * np.sqrt(term) def cv_simulation(V_g, polarization_func, direction='up'): """ C-V特性仿真 基于文件公式(4.1)-(4.8) """ C_total = np.zeros_like(V_g) psi_s0, E_f0 = 0.1, 0.1 # 初始猜测值 for i, V_g_val in enumerate(V_g): def equations(vars): psi_s, E_f = vars # 公式(4.3): P(E_f) = -Q_s(ψ_s) P = polarization_func(E_f, direction) Q_s_val = Q_s(psi_s) eq1 = P + Q_s_val # 公式(4.2): V_F = E_f * t_F V_f = E_f * params['t_F'] # 公式(4.3): V_I = -Q_s(ψ_s) * EOT / (ε_0 ε_i) V_i = -Q_s_val * params['EOT'] / (epsilon_0 * params['epsilon_i']) # 公式(4.1): V_G = ψ_s + V_I + V_F eq2 = V_g_val - (psi_s + V_i + V_f) return [eq1, eq2] # ============================================================= # 改进的牛顿法解器 # ============================================================= def newton_solver(equations, initial_guess, max_iter=50, tol=1e-6): """ 牛顿法解非线性方程组 equations: 函数,返回方程残差 [f1, f2] initial_guess: 初始猜测 [x0, y0] 返回: 解 [x, y] """ x, y = initial_guess for i in range(max_iter): # 计算函数值 f = np.array(equations([x, y])) # 检查收敛 if np.linalg.norm(f) < tol: return x, y # 数值微分计算雅可比矩阵 dx = 1e-6 J = np.zeros((2, 2)) # df1/dx f1_x = equations([x + dx, y])[0] f1 = equations([x, y])[0] J[0, 0] = (f1_x - f1) / dx # df1/dy f1_y = equations([x, y + dx])[0] J[0, 1] = (f1_y - f1) / dx # df2/dx f2_x = equations([x + dx, y])[1] f2 = equations([x, y])[1] J[1, 0] = (f2_x - f2) / dx # df2/dy f2_y = equations([x, y + dx])[1] J[1, 1] = (f2_y - f2) / dx # 解线性方程组 J * delta = -f try: delta = np.linalg.solve(J, -f) except np.linalg.LinAlgError: # 奇异矩阵,使用伪逆 delta = np.linalg.pinv(J) @ (-f) # 更新解 x += delta[0] y += delta[1] # 应用物理约束 x = max(min(x, 20.0), -20.0) # 表面势约束 y = max(min(y, 1e7), -1e7) # 铁电电场约束 # print(i,np.linalg.norm(f)) # 未收敛 raise RuntimeError(f"牛顿法未在 {max_iter} 次迭代内收敛") # ============================================================= # 改进的 C-V 特性仿真 # ============================================================= def cv_simulation(V_g, polarization_func, direction='up'): C_total = np.zeros_like(V_g) # 存储前一个解用于初始猜测 prev_psi_s, prev_E_f = 0.1, 0.1 for i, V_g_val in enumerate(V_g): def equations(vars): psi_s, E_f = vars P = polarization_func(E_f, direction) Q_s_val = Q_s(psi_s) V_f = E_f * params['t_F'] V_i = -Q_s_val * params['EOT'] / (epsilon_0 * params['epsilon_i']) eq1 = P + Q_s_val eq2 = V_g_val - (psi_s + V_i + V_f) return [eq1, eq2] try: # 基于栅压的智能初始猜测 if i == 0: # 第一点:据栅压符号猜测 if V_g_val > 0: initial_guess = [0.5, 1] # 正栅压时表面势和电场较大 else: initial_guess = [-0.5, 1] # 负栅压时表面势和电场较小 else: # 使用前一个解作为初始猜测 initial_guess = [prev_psi_s, prev_E_f] # 使用牛顿法解 psi_s, E_f = newton_solver(equations, initial_guess) prev_psi_s, prev_E_f = psi_s, E_f # 使用公式(4.5)计算半导体电容 term = (params['n_i']**2 / params['N_D']**2) * (-np.exp(-beta * psi_s) + 1) + \ (np.exp(beta * psi_s) - 1) # 添加数值保护 if term < 0: term = 1e-6 # 避免负值 C_si = (params['A_f'] * epsilon_0 * params['epsilon_i']) / (np.sqrt(2) * params['L_D']) * np.sqrt(term) # 铁电层电容 C_f = params['A_f'] * epsilon_0 * params['epsilon_F'] / params['t_F'] # 绝缘层电容 C_i = params['C_i'] # 计算总电容 C_total[i] = 1 / (1/C_f + 1/C_i + 1/C_si) except Exception as e: if i > 0: C_total[i] = C_total[i-1] else: C_total[i] = 1e-12 print(f"CV解失败: V_g={V_g_val}, 错误: {e}") return C_total # ============================================================= # 改进的 I-V 特性仿真 # ============================================================= def solve_surface_potential(V_g, V, polarization_func, direction): def equations(vars): psi_s, E_f = vars P = polarization_func(E_f, direction) Q_s_val = Q_s(psi_s) V_i = -Q_s_val * params['EOT'] / (epsilon_0 * params['epsilon_i']) V_f = E_f * params['t_F'] eq1 = P + Q_s_val eq2 = V_g - (psi_s + V_i + V_f) return [eq1, eq2] # 基于栅压的智能初始猜测 if V_g > 0: initial_guess = [0.5, 1] # 正栅压 else: initial_guess = [-0.5, -1] # 负栅压 try: psi_s, E_f = newton_solver(equations, initial_guess) return psi_s except: # 回退到物理合理的默认值 return abs(params['psi_b']) def pao_sah_integrand(psi, V): """ 改进的Pao-Sah被积函数,添加数值保护 """ # 计算公式(21)中的内部表达式 term1 = np.exp(beta * psi) - beta * psi - 1 term2 = (params['n_i']**2 / params['N_D']**2) * np.exp(beta * V) * \ (np.exp(-beta * psi) + beta * psi * np.exp(-beta * V) - 1) # 添加数值保护 total_term = term1 + term2 # print(total_term) if total_term <= 0: # 物理上应为正,处理数值误差 total_term = 1e-10 # 计算电场 ξ(ψ,V) prefactor = (2 * params['N_D'] * k * T) / (epsilon_0 * params['epsilon_Si']) xi = np.sqrt(prefactor * total_term) # 公式(20)被积函数 numerator = (params['n_i']**2 / params['N_D']) * np.exp(-beta * (psi - V)) # 避免除零 if xi < 1e-10: return 0 return numerator / xi def pao_sah_integral(V_ds, V_g, polarization_func, direction): """ 改进的Pao-Sah双积分,使用自适应积分 """ # 使用自适应步长 num_V_points = max(10, int(V_ds * 10)) # 步长随V_ds自适应 V_points = np.linspace(0, V_ds, num_V_points) integral_val = 0 for i in range(len(V_points)): V_val = V_points[i] try: psi_s = solve_surface_potential(V_g, V_val, polarization_func, direction) except: psi_s = abs(params['psi_b']) if psi_s > params['psi_b']: # 自适应步长:在变化大的区域增加点数 num_psi_points = max(20, int((psi_s - params['psi_b']) * 100)) psi_points = np.linspace(params['psi_b'], psi_s, num_psi_points) #print(psi_s) # 计算被积函数值 integrand_vals = [] for psi in psi_points: val = pao_sah_integrand(psi, V_val) # 处理异常值 if not np.isfinite(val) or val < 0: val = 0 integrand_vals.append(val) inner_integral = trapz(integrand_vals, psi_points) # 计算dV (梯形法) if i == 0: dV = V_points[1] - V_points[0] elif i == len(V_points) - 1: dV = V_points[-1] - V_points[-2] else: dV = (V_points[i+1] - V_points[i-1]) / 2 integral_val += inner_integral * dV return integral_val def iv_simulation(V_g, V_ds, polarization_func, direction='up'): """ I-V特性仿真 基于文件公式(4.9): I_D = q μ (W/L) * [Pao-Sah双积分] """ I_d = np.zeros_like(V_ds) for i, V_ds_val in enumerate(V_ds): try: integral = pao_sah_integral(V_ds_val, V_g, polarization_func, direction) I_d[i] = q * params['mu'] * (params['W'] / params['L']) * integral except Exception as e: if i > 0: I_d[i] = I_d[i-1] else: I_d[i] = 1e-12 print(f"IV解失败: V_ds={V_ds_val}, 错误: {e}") return I_d # ============================================================= # 主程序 # ============================================================= if __name__ == "__main__": # 设置扫描电压范围 V_g_cv = np.linspace(-1.6, 0.4, 30) # C-V栅压扫描 V_ds_iv = np.linspace(-0.5, 0, 20) # I-V漏压扫描 V_g_iv = 3.0 # I-V固定栅压 V_ds_transfer = 0.1 # 转移特性固定漏压 print("开始C-V仿真 (上扫)...") C_up = cv_simulation(V_g_cv, external_polarization, direction='up') print("开始C-V仿真 (下扫)...") C_down = cv_simulation(V_g_cv[::-1], external_polarization, direction='down') C_down = C_down[::-1] # 反转结果以匹配电压顺序 print("开始I-V仿真 (ON状态)...") I_d_on = iv_simulation(V_g_iv, V_ds_iv, external_polarization, direction='up') print("开始I-V仿真 (OFF状态)...") I_d_off = iv_simulation(V_g_iv, V_ds_iv, external_polarization, direction='down') print("开始转移特性仿真...") V_g_transfer = np.linspace(0, 5, 20) I_d_transfer = np.zeros_like(V_g_transfer) for i, V_g_val in enumerate(V_g_transfer): I_d_transfer[i] = iv_simulation(V_g_val, [V_ds_transfer], external_polarization, direction='up')[0] # 绘制结果 plt.figure(figsize=(15, 12)) # 1. P-E曲线(新增子图) plt.subplot(2, 3, 1) plt.plot(E_up, P_up, 'b-', label='Up Sweep') # 直接使用原始数据 plt.plot(E_down, P_down, 'r--', label='Down Sweep') plt.xlabel('Electric Field (MV/m)') # 单位已正确 plt.ylabel('Polarization (C/m²)') # 单位已正确 plt.title('Ferroelectric P-E Curve') plt.grid(True) plt.legend() # 2. C-V曲线(原1号子图调整到2号位置) plt.subplot(2, 3, 2) plt.plot(V_g_cv, C_up * 1e6, 'b-',label='Up Sweep') plt.plot(V_g_cv, C_down * 1e6, 'r--', label='Down Sweep') plt.xlabel('Gate Voltage (V)') plt.ylabel('Capacitance (μF)') plt.title('MFIS Capacitor C-V Characteristics') plt.grid(True) plt.legend() # 3. I-V曲线 (ON/OFF状态)(原2号子图调整到3号位置) plt.subplot(2, 3, 3) plt.plot(V_ds_iv, I_d_on * 1e6, 'g-', label='ON State') plt.plot(V_ds_iv, I_d_off * 1e6, 'k--', label='OFF State') plt.xlabel('Drain Voltage (V)') plt.ylabel('Drain Current (μA)') plt.title(f'FeMFET I-V Characteristics (V_g = {V_g_iv}V)') plt.grid(True) plt.legend() # 4. 不同栅压下的I-V曲线(原3号子图调整到4号位置) plt.subplot(2, 3, 4) V_g_list = [1.0, 2.0, 3.0, 4.0] for V_g_val in V_g_list: I_d = iv_simulation(V_g_val, V_ds_iv, external_polarization, direction='up') plt.plot(V_ds_iv, I_d * 1e6, label=f'V_g = {V_g_val}V') plt.xlabel('Drain Voltage (V)') plt.ylabel('Drain Current (μA)') plt.title('FeMFET I-V at Different Gate Voltages') plt.grid(True) plt.legend() # 5. 转移特性曲线(原4号子图调整到5号位置) #plt.subplot(2, 3, 5) #plt.semilogy(V_g_transfer, np.abs(I_d_transfer) * 1e6, 'mo-') #plt.xlabel('Gate Voltage (V)') #plt.ylabel('Drain Current (μA)') #plt.title(f'Transfer Characteristics (V_ds = {V_ds_transfer}V)') #plt.grid(True) plt.tight_layout() # 保存图像 output_dir = r"D:\code\phi" if not os.path.exists(output_dir): os.makedirs(output_dir) plt.savefig(os.path.join(output_dir, "mfisfet_simulation_results.png"), dpi=300) print(f"结果已保存至: {os.path.join(output_dir, 'mfisfet_simulation_results.png')}") plt.show() # 保存数据 np.savetxt(os.path.join(output_dir, "cv_data.txt"), np.column_stack((V_g_cv, C_up * 1e6, C_down * 1e6)), header="V_g(V) C_up(μF) C_down(μF)") np.savetxt(os.path.join(output_dir, "iv_data.txt"), np.column_stack((V_ds_iv, I_d_on * 1e6, I_d_off * 1e6)), header="V_ds(V) I_on(μA) I_off(μA)") print("仿真数据已保存")现在解器有巨大问题,积分和解方程都不收敛。能不能用matlab写一个。matlab有没有能直接解方程的工具
06-14
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值