# Numerical solution for multi-step PDEs

I am trying to numerically solve the following PDEs in MMA (v12.0) $$\partial_u f(x,u) = \left\{\begin{matrix}&-a f(x,u) + b \, \partial_x^2 f(x,u), & 0

with free boundary conditions $$f(x\to\pm\infty, u)=0$$. The inputs are the initial condition $$f(x,0)$$, the number of "rounds" $$n$$ (each round consists of two steps), parameters $$a$$ and $$b$$, together with the function $$r$$ (a specific example is $$r(x;\alpha) = -\alpha x^2$$) and a set of parameters $$\alpha_1,\alpha_2,\cdots,\alpha_n$$. The number of rounds is going to be around $$n \lesssim 20$$. The output $$f_{\rm sol}$$ ideally gives the value for all $$0 < u < n(T+\tau)$$.

The problem is that I don't know how to automatize the process of evaluating the solution given a value for $$n$$, without having to explicitly write a code to solve for each step in each round, or having to explicitly enter the equations repeatedly (since they have the same form and only different parameters); the explicit solution can of course be done using NDSolve, for example. Any input regarding this is appreciated (even if it works for a fixed $$n$$, for example $$n=10$$).

My (half-)Attempt

Since these equations are simple (heat equation for the first steps, linear growth for the second steps), one can write explicit expressions for their solutions in each step, and then pass the value of the function at the end of the interval into the solution of the next step.

For example, for the first line, the solution is given as the convolution of the initial condition with the heat kernel as $$f(x, u) = e^{-au} \int dx_0 \, \underbrace{\frac{\exp\left[ -\frac{(x-x_0)^2}{4bu}\right]}{\sqrt{4\pi b u}}}_{H(x-x_0,u)} \, f(x_0,0), \qquad \text{for} \qquad 0 and for the second step, the solution is the exponential growth given by $$‌ f(x, u) = \exp\left[ r(x;\alpha_1) (u-T) \right] \, f(x,T) \qquad \text{for} \qquad T

So what I'm thinking is to define an MMA function for the first step as _roughly, this exact syntax doesn't work - I don't know how to implement this procedure _

H[x_,u_,b_] := Exp[-x^2/(4 b u)]/Sqrt[4 Pi b u];
FirstStep[g_[x_],x_,u_,a_,b_] := Exp[-a T] Convolve[g[x0], H[x0,u,b], x0, x];
SecondStep[g_[x_],x_,u_,r_[x_,alpha_]] := Exp[r[x,alpha]*(u-u_0)] g[x];
OneRound[g_[x_],x_,u_,a_,b_,r_[x_,alpha_]]:= SecondStep[FirstStep[g,x,u,a,b],x,u,r]


OneRound is the consequent application of SecondStep after FirstStep; I need to apply OneRound to the initial condition, say given by g[x0_]:= PieceWise[ {-x0^2+4, -2<x0<2}, 0].

I vaguely think one could use Fold (instead of Nest, such that the arguments $$\alpha_1,\alpha_2,\cdots$$ can be passed to the function) to apply FirstStep and SecondStep to f[x,0] a number of $$n$$ times, maybe something like

This approach, however, would only give the final form $$f(x,n(T+\tau))$$ and not the solution at all times - but it'd still be useful. Any suggestion for implementing this approach would be very helpful, too.

• If I understand your question correctly, I think you just need a Piecewise? May 18, 2022 at 1:28
• @xzczd could you please show me how? I need to solve the pde for each interval, and then use the value of the solution at the end of that interval as the initial value for the pde in the next interval, etc. May 18, 2022 at 1:40
• Just search "NDSolve Piecewise" in this site, you'll see tons of examples. If you still feel confused, please add a specific example to the question so we can play with it. May 18, 2022 at 1:51
• @xzczd it seems that with Piecewise I'd have to repeatedly enter the rhs for every time interval - which for large number of such intervals would be cumbersome. May 18, 2022 at 3:06
• Just use a Table. May 18, 2022 at 4:40