| 查看: 563 | 回复: 0 | ||
[求助]
Matlab解ODE方程组 出现singular
|
|
Hi, 我需要求解一个ode方程组,目前用的matlab的ode23s。但是求解不出来,算出第一个数值后,解不了第二步,显示错误信息 "Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552643e-19) at time t." 或者 "Matrix is close to singular or badly scaled. Results may be inaccurate. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552643e-19) at time t." 能否麻烦大家帮我看一下问题在哪里?非常极其感谢! Matlab程序如下: global F R T e0 Z_Na Z_Cl Z_OH Z_H C_Na_A C_Cl_A C_OH_A C_H_A Cond_A i U_A E_A FixedCharge LA D_H D_OH e_r epsilon C_Na_A2 C_Cl_A2 C_H_A2 C_OH_A2 F=96485; R=8.3144621; T=25+273.15; e0=8.85419/10^12; Z_Na=1; Z_Cl=-1; Z_OH=-1; Z_H=1; %%% solution = 1 M NaCl + 0.1 M NaOH in the BULK C_Na_A=1000+100; C_Cl_A=1000+1e-10; C_OH_A=100; C_H_A=1e-10; Cond_A=10.1; i=1000; U0=0; E_A=i/Cond_A; FixedCharge=2*1000; LA=1E-4; D_H=5.94/10^10; D_OH=3.47/10^10; e_r=20; epsilon=e_r*e0; U_A=-R*T/F*asinh(-FixedCharge/2/(C_Na_A+C_H_A))+U0; C_Na_A2=C_Na_A*exp(-Z_Na*F/R/T*U_A); C_Cl_A2=C_Cl_A*exp(-Z_Cl*F/R/T*U_A); C_OH_A2=C_OH_A*exp(-Z_OH*F/R/T*U_A); C_H_A2=C_H_A*exp(-Z_H*F/R/T*U_A); %%%%%%%%%%%%%%%%%%%% y0=[E_A;U_A;C_OH_A2]; [p Y]=ode23s(@odefun,linspace(-LA,0,10000),y0); %%%%%%%%%%%%%%%%%%%% 方程组function如下: function [dydx]= odefun(x,y) global F R T Z_Na Z_Cl Z_OH Z_H C_Na_A2 C_Cl_A2 i U_A FixedCharge D_OH epsilon D_H C_Na C_Cl % y(1)=Electric field, y(2)=Electric potential, y(3)=OH- concentration dydx=zeros(3,1); C_Na=C_Na_A2*exp(-Z_Na*F/R/T.*(y(2)-U_A)); C_Cl=C_Cl_A2*exp(-Z_Cl*F/R/T.*(y(2)-U_A)); dydx(1)=(F/epsilon)*(C_Na-C_Cl+10^(-8)./y(3)-y(3)+FixedCharge); dydx(2)=-y(1); J_H=-D_H*(-10^(-8)./(y(3)^2)*dydx(3)+Z_H*F/R/T*10^(-8)./y(3).*dydx(2)); dydx(3)=-((i/F-Z_H*J_H)/Z_OH/D_OH)-Z_OH*F/R/T.*y(3).*dydx(2); end |
» 猜你喜欢
焦虑
已经有13人回复
一志愿985,本科211,0817化学工程与技术319求调剂
已经有8人回复
0817 化学工程 299分求调剂 有科研经历 有二区文章
已经有9人回复
346求调剂[0856]
已经有5人回复
285化工学硕求调剂(081700)
已经有11人回复
267一志愿南京工业大学0817化工求调剂
已经有9人回复
274求调剂
已经有6人回复
求材料调剂
已经有10人回复
本科郑州大学物理学院,一志愿华科070200学硕,346求调剂
已经有4人回复
国自科面上基金字体
已经有8人回复













回复此楼