New answers tagged finite-difference
1
Verification of order of convergence of Implicit Euler Method to numerically solve Black-Scholes PDE
You forgot to multiply the norm of the difference between the numerical solution and the exact one by the discretization step. In your case it is enough to divide the err_n by Nt when computing the order of time discretization, and by NS when computing the order of space discretization. I got a perfect first order accuracy for the time discretization, in the ...
0
So if I get you right, what you are recommending is to reorganise the system of equations like this:
$$
\displaystyle \frac{\partial c}{\partial t } = \frac{-D_{e}}{h^2} \begin{bmatrix}
0 & 0 & \dots & 0 & 0 \\
1 & -2 & 1 & \dots ...
3
You did nothing wrong, the numbers just are what they are. Plotting the error for the backwards difference gives the plot
which is as expected for a method of order one with step size about $h=0.01$. For the second order methods one expects an error of the magnitude $h^2=10^{-4}$, the plot below confirms this
2
Yes. That's all there is to the stability condition.
Taking the material properties - shear modulus ($\mu$), bulk modulus ($\kappa$) and density ($\rho$) - into account, the global critical time step is evaluated as the minimum of the critical time step for each element ($\Delta t^e$)
$\Delta t^e = CFL * h^e / c_{\kappa}$
where CFL is the Courant-Friedrichs-...
2
Let's reconstruct this from first principles:
Defining the ODE system
In the method-of-lines discretization you solve an ODE system $\dot U=F(U)$, $U=(U_0,U_1,...,U_{M+1})$, $U_k(t)=u(x_k,t)$, and similarly $F=(F_0,F_1...,F_{M+1})$. Because of the boundary conditions
$$
u(x, 0) = 40 · x^2 · (1 - x) / 3
\\
u(0, t) = u(1, t) = 0
$$
$U_0=U_{M+1}=0$ and ...
2
There is really no such thing as a good finite difference equivalent to an operator. In the earliest days of scientific computing, the thought was that each differential operator would be replaced by some finite difference expression, and that finite difference operator would be the most accurate one available, usually a central difference.
I believe that ...
0
Some remarks about the code, perhaps cleaning these up already solves the question, especially the initialization error:
Floating point errors are a thing, so if you want to get the expected results, give some wiggle room (0.01 is arbitrary, I would not expect defects above 1e-12 in this place)
N = int((x1 - x0) / h + 0.01)
M = int((t1 - t0) / tau + 0.01)...
Top 50 recent answers are included
