enter image description hereI want to solve for s[ i ] ,L[ i ] using the two equations f1 and f2, i know the values of s, L at time t=0, and i want s,L values at different time steps. My second equation involves s, L value from current time step [i] and from previous time step [i-1].I tried to code this but couldn' t get, please help me to solve this (Equations picture is attached)
from scipy.optimize import fsolve
import numpy as np
import math as m
q0=0.0011
thetas,thetai,thetar=0.43,0.1,0.05
ks=0.0022#m/hr
psib=-0.15# m
lamda=1
eta=2+3*lamda
ki=8.676*10**(-8)
si=0.13157
tmin = 0
tmax = 30
dt = 0.01
t = np.arange(tmin,tmax+dt,dt)
n=np.size(t)
S=0.02/24
s=np.zeros(n)
L=np.zeros(n)
# s[0]=0.13157
# L[0]=0
delta=-0.1001
def equations(p):
s, L = p
f1=(ks*s[i]**(3+(2/lamda))-(psib/(1-eta))*(((ki*si**(-1/lamda))-(ks*s[i]**(3+(1/lamda))))/L[i])-q0)
f2=(L[i]*(s[i]*(thetas-thetar))+S*t[i]*dt*0.5*((L[i]*m.exp(-delta*psib*(-1+s[i]**(-1/lamda))))+(L[i-1]*m.exp(-delta*psib*(-1+s[i-1]**(-1/lamda)))))-(q0-ki)*t[i])
return(f1,f2)
for i in range(n):
if i==0:
s[i]=si
L[i]=0
s,L=fsolve(equations,([0.19,0.1]))
theta=s*(thetas-thetar)+thetar
print("For soil 2 @ t={} soil moisture is {} & wetting front depth in cm is {}".format(t[i],theta,L*100))
else:
s,L=fsolve(equations,([0.19,0.1]))
theta=s*(thetas-thetar)+thetar
print("For soil 2 @ t={} soil moisture is {} & wetting front depth in cm is {}".format(t[i],theta,L*100))
