Stack Exchange Network

Stack Exchange network consists of 175 Q&A communities including Stack Overflow, the largest, most trusted online community for developers to learn, share their knowledge, and build their careers.

Visit Stack Exchange
5

The upwind scheme is fine. The code below is improved. There were various syntax errors on data types (scalars, vectors, matrices). %clear workspaces clear clc % define variables xmin=-10; % minimum value of x xmax=10; % maximum value of x N=100; % no. nodes - 1 dt=0.005; % timestep t=0; % time tmax=10; % maximum ...


4

There are several ways to think about approximating $\tan x$ as a rational function. Taking the log-derivative of $\sin x=x\prod_{k\ge 1}\left(1-\frac{x^2}{k^2\pi^2}\right)$ gives $$\cot x=\frac{1}{x}-2\sum_{k\ge 1}\frac{x}{k^2\pi^2-x^2}.$$Keeping only the $k=1$ term,$$\tan x\approx\frac{x(\pi^2-x^2)}{\pi^2-3x^2}.$$This is a reasonable approximation if $|x|\...


4

You can also directly apply Dirichlet's test (the integral of the numerator is zero when taken over a period, therefore uniformly bounded on all finite intervals). We have $$\cos^{2 k} u = 2^{-2 k} \binom {2 k} k + 2^{1 - 2 k} \sum_{1 \leq j \leq k} \binom {2 k} {k - j} \cos 2 j u, \\ \sum_{k \geq 0} \frac {(i \alpha \cos u)^{2 k}} {(2 k)!} = J_0(\alpha) + ...


4

I don't see an elegant way to do this so we will apply brute force starting from this formula $$w_k=\frac{a_n}{a_{n-1}}\frac{\int_a^bw(x)p_{n-1 }(x)^2dx}{p_n^{\prime}(x_k)p_{n-1}(x_k)}$$ In our case we write it as $$v_k=\frac{b_m}{b_{m-1}}\frac{\int_0^1\sqrt y\,q_{m-1}(y)^2dy}{q_m^{\prime}(y_k)q_{m-1}(y_k)}$$ Now $b_m$ is the leading coefficient of $q_m(y)=\...


4

Hints: (1) $f$ extends to a continuous function $[0,\infty)\to\mathbb{R}$. (2) $\lim_{x\to\infty}f(x)=0$.


3

You know by Taylor expansions that $u(t+h)=u(t)+hu'(t+h/2)+\frac{h^3}6u'''(t+h/2)+O(h^5)$. Insert the ODE for the first derivative and compare with the numerical method \begin{align} u(t+h)&=u(t)+hf(t+h/2, u(t+h/2))+\frac{h^3}6u'''(t+h/2)+O(h^5)\\ U_{+1}&=U+hf(t+h/2, \tilde U)\\[1em]\hline u(t+h)-U_{+1}&=u(t)-U+h\left[f(t+h/2, u(t+h/2))-f(t+h/2, \...


3

A defining property of the outer product is its anti-symmetry $$ a∧b = -b∧a, $$ which implies that $$ 2\,a∧a=0. $$


3

Let us diagonalize the matrix $A$ such that ${\bf u}_t + A {\bf u}_x = {\bf 0}$ as follows: \begin{aligned} A &= \begin{pmatrix}0 & a\\ b & 0\end{pmatrix} \\ &= \begin{pmatrix}-\sqrt{a/b} & \sqrt{a/b}\\ 1 & 1\end{pmatrix} \begin{pmatrix}-\sqrt{a/b} & 0\\ 0 & \sqrt{a/b}\end{pmatrix} \frac12 \begin{pmatrix}-\sqrt{b/a} & 1\\ \...


3

Your approach works, because it yields another ODE: Let $z(t)=x(T-t)$. Then $$ \frac{\partial }{\partial t}z(t)=\frac{\partial }{\partial t}x(T-t)\\ =-x'(T-t)=-f(T-t,x(T-t))\\ =-f(T-t,z(t)). $$


3

In fact, that matrix does not have a dominant eigenvalue. You have that $|\lambda_1| = |\lambda_2|=|\lambda_3| =1$ and so there is no reason why the method should converge, unless you take as initial approximation an eigenvector associated with $\lambda_1$. The exact values for $\lambda_2, \lambda_3$ are $-\frac 12 \pm \frac{\sqrt{3}}{2}$.


3

You should check out chapter 3 in Demmel's text. Let me summarize some of the results for dense least-squares, for $A$ an $m\times n$ matrix: $A$ is well-conditioned: In this case, using the normal equations is around as accurate as other methods and is also the fastest, so using them is fine, even though they are not numerically stable. $A$ is not well-...


2

Start with the Taylor expansion $$ u(x,t+k)=u(x,t)+u_t(x,t)k+\frac12u_{tt}(x,t)k^2+... $$ Use the PDE to convert time into space derivatives, $u_t=-[f(u)]_x$, then \begin{align} u_{tt}&=-[f(u)]_{xt}=-[f'(u)u_t]_x\\ &=[f'(u)[f(u)]_x]_x\\ \end{align} Now realize this derivative structure using symmetric divided differences and a half-step scheme. Let ...


2

You get an automatic Newton solver by using the BVP solver scipy.integrate.solve_bvp. You program your equation as usual def eqn(x,u): return [ u[1], (1+u[1]**2)**1.5*(T*u[0]+0.5*w*x*(x-L))/K ] def bc(u0,uL): return [ u0[0], uL[0] ] x = np.linspace(0,L,21); a = 0.1 u = [ a*x*(L-x), a*(L-2*x) ] res = solve_bvp(eqn, bc, x, u, tol = 1e-12); # the solution is ...


2

The interval is quite simply a consequence of following the standard proof of Picard-Lindelöf. As the Lipschitz constant is globally $L=1$, one does not need a restriction in the $y$ direction. In the next step, the Picard iteration is considered on $C([−ϵ,ϵ])$ where it has a Lipschitz constant as a mapping on a function space of $Lϵ=ϵ$, $$ \bigl|P[y_1](t)...


2

The value of $y'(a)$ depends on the method you use and the step size. However, because your system is linear, solving for $y'(a)$ is just a single step linear interpolation/extrapolation. Alternatively, you can keep track of $y'(a)$ during your iterations (since the DE is linear, every step of the approximate solution is going to look like $y=y_0+y'(a)y_1$)...


2

This is not an identity, but just an approximation for $\tan$. If you look at their graphs, you'll notice that the two functions line up very, very well though for $|x|\leq\pi$.


2

Yes, in $c_1h^p+c_2h^{p+1}+...$, the second term will dominate the first one for $h>\frac{c_1}{c_2}$. In numerical applications, the many steps required by smaller step sizes eventually accumulate floating point noise sufficient to dominate the truncation error, so that a loglog plot of error vs. step size has a V shape with a fuzzy left leg, a middle ...


2

Let me reformulate your assertion under the precise form : $$\lim_{x\to \infty}\left(\sqrt[n]{x^n+q(x)}-x\right)=\dfrac{1}{n}c_{n-1} \ \ \ \text{where} \ \ \ q(x):=c_{n-1}x^{n-1}+\cdots+c_0 \tag{1}$$ and prove it. Indeed : $$\sqrt[n]{x^n+q(x)}=\sqrt[n]{x^n\left(1+\dfrac{q(x)}{x^n}\right)}=x\sqrt[n]{1+\dfrac{q(x)}{x^n}}.\tag{2}$$ Now, let us use the first ...


2

It is easy to check that you get from the Euler method to a second order method by taking the slope at the midpoint as in the explicit midpoint or modified Euler method $$y(x+Δx)=y(x)+f(x+\tfrac12Δx, y(x)+\tfrac12f(x,y(x))Δx)Δx$$ or by combining two slopes as in the explicit trapezoidal or improved Euler or Heun's 2nd order method $$y(x+Δx)=y(x)+\tfrac12f(x,...


2

Calculate backwards, with partial integration: \begin{align} \int_a^b(x-a)(x-b)f''(x)dx &=[(x-a)(x-b)f'(x)]_a^b-\int_a^b(2x-a-b)f'(x)dx\\ &=0-0-[(2x-a-b)f(x)]_a^b+\int_a^b2f(x)dx\\ &=-(b-a)(f(b)+f(a))+2\int_a^bf(x)dx \end{align} so that indeed the claimed identity holds.


2

No, it is not. Take $f(x)=1$ for all $x$, and take $g(x)=1$ if $x$ is rational and $g(x)=-1$ if $x$ is irrational. $f$ and $f\circ g$ are both constant, but $g$ is nowhere continuous.


2

Assuming IEEE754 double representation (so there is a sign bit on exponent and mantissa), the largest $P$ is $\approx 2^{2^{11-1}}=2^{1024}$, so $x^{64}\approx 2^{1024}$ and hence $x\approx 2^{16}=65536$. This ignores all subtleties about precision of the result depending on the machine implementation, order of operations, etc.


2

As you want to compute the inverse, let us start from that end: Consider $$g(x)=\begin{cases}3x^5-10x^3+30x&-1\le x\le 1\\ 15x+8&x\ge1\\ 15x-8&x\le -1\end{cases} $$ Then $g$ is $C^2$ and has no critical points, hence has a $C^2$ inverse $f$ that looks like a line plus a sigmoidal. To match your demands for $x\ll 0$ and $x\gg 0$, we consider a ...


2

To solve 2D quasilinear systems of conservation laws $$ {\bf u}_t + {\bf A}({\bf u})\, {\bf u}_x + {\bf B}({\bf u})\, {\bf u}_y = {\bf 0} $$ numerically, various strategies can be followed: Implement a 2D finite-volume scheme, such as the 2D Lax-Friedrichs method \begin{aligned} {\bf u}_{i,j}^{n+1} &= \frac{{\bf u}_{i-1,j}^{n} + {\bf u}_{i+1,j}^{n} + {\...


2

I think that the process you're envisioning, if you did it rigorously, becomes equivalent to Picard's iterative process: Given an initial value problem $y'(t)= f(t,y)$, $y(0) = y_0$, the sequence of functions defined by \begin{align*} \phi_0(t) &= y_0 \\ \phi_{k+1}(t) &= y_0 + \int_0^t f(s,\phi_k(s))\,ds \end{align*} converges to a solution $\...


2

In the last line of your code, you have h/2. It should be h/3. You also are using the trapezoid weights instead of the simpson's weights. In fact, I can't figure out why your two results are different at all, since the calculations in the last two lines are identical.


2

Another point of view is the sampling theorem, as the integrated function is periodic and integrated over 2 periods. The limit frequency of $\sin^2x =\frac12(1-\cos2x)$ is $2$, so with 4 sub-intervals you are at the minimal sampling frequency. If you write $S(h)=\frac{4T(h)-T(2h)}3$ as per Richardson extrapolation, then the term $T(2h)$ is under-sampled with ...


2

You are missing an important part of the story... After you multiply the equation by a smooth test function, you extend it by density to some larger (Sobolev) space. The finite element spaces will be subspaces of that Sobolev space, not of $C^{\infty}_c(\Omega)$, so they are not supposed to be differentiable across elements.


1

No! Note the bar on $\bar{\mathcal{F}}$. You rewrite every appearance of $\Delta_\delta y(x)=\frac{y(x+\delta)-y(x)}{\delta}$, $\Delta_\delta^2 y(x)=\frac{y(x+2\delta)-2y(x+\delta)+y(x)}{\delta^2}$, etc., so an equation of the form $$ \mathcal{F}(y(x),\Delta_\delta y(x), \Delta_\delta^2 y(x),\dots;x)=0 $$ becomes an equation involving only $x, y(x), y(x+\...


1

The best course of action is to reduce to a system of first order differential equations and use a standard ode solver on that system. $$ \begin{cases} y' = z\\ z' = -\frac{k}{m} (y-x) +\frac fm \\ x' = \frac kb (y-x) \end{cases} $$ Edit: With the set of parameters given in the comments below, Wolfram returns a numerical solution with a maximum error in ...


Only top voted, non community-wiki answers of a minimum length are eligible