6
$\begingroup$

I don't quite understand the process of solving differential equations by MATLAB. It seems that it doesn't need the explicit function to specify the required solution, but only needs to input the vector formed by its derivative into ode45 to solve it. I want to implement the same code in Mathematica, but I met with difficulties.

This is the MATLAB code I want to convert to Mathematica code:

clear all; close all; 
L=80; N=256; 
x=L/N*[-N/2:N/2-1]; 
k=(2*pi/L)*[0:N/2-1 -N/2:-1].'; 
%init conditons 
u=2*sech(x); ut=fft(u); 
vt=zeros(1,N); uvt=[ut vt]; 
%Solve
a=1; t=0:0.5:20; 
[t,uvtsol]=ode45('wave1D',t,uvt,[],N,k,a); 
usol=ifft(uvtsol(:,1:N),[],2); 
%Plot 
p=[1 11 21 41]; 
for n=1:4 
subplot(5,2,n) 
plot(x,usol(p(n),:),'k','LineWidth',1.5), xlabel x, ylabel u 
title(['t=' num2str(t(p(n)))]), axis([-L/2 L/2 0 2]) 
end 
subplot(5,2,5:10) 
waterfall(x,t,usol), view(10,45) 
xlabel x, ylabel t, zlabel u, axis([-L/2 L/2 0 t(end) 0 2]) 

File wave1D.m code:

function duvt=wave1D(t,uvt,dummy,N,k,a) 
ut=uvt(1:N); vt=uvt(N+[1:N]); 
duvt=[vt; -a^2*(k).^2.*ut]; 

This is what I've tried:

ClearAll["`" ~~ _]
L = 80;
N0 = 256;
x = L/N0*Range[-N0/2, N0/2 - 1];
k = (2*Pi/L)*Join[Range[0, N0/2 - 1], Range[-N0/2, -1]];

u = 2*Sech[x];
ut = Fourier[u, FourierParameters -> {-1, 1}];
vt = ConstantArray[0, N0];
uvt = Join[ut, vt];

a = 1;
ufl = Table[uf[i][0], {i, 1, N0}];
vfl = Table[vf[i][0], {i, 1, N0}];
dufl = Table[D[uf[i][t], t], {i, 1, N0}];
dvfl = Table[D[vf[i][t], t], {i, 1, N0}];
eqs1 = Thread[Equal[dufl, vt]];
eqs2 = Thread[Equal[dvfl, -a^2 k^2 ut]];
int1 = Thread[Equal[ufl, ut]];
int2 = Thread[Equal[vfl, vt]];
rs = NDSolveValue[{eqs1, eqs2, int1, int2}, {ufl, vfl}, {t, 0, 20}];

But I can't get the correct result, I don't know how to fix it. Please help!

$\endgroup$
10
  • 1
    $\begingroup$ [t,uvtsol]=ode45('wave1D',t,uvt,[],N,k,a); Error using feval Unrecognized function or variable 'wave1D'. what is wave1D your code does not run for me on Matlab. Do you have some special toolbox installed that defines this? Any way, asking folks here to translate Matlab code to Mathematica is probably not covered by the mission of this forum, as this forum is really for Mathematica programming questions. $\endgroup$
    – Nasser
    Dec 4, 2023 at 1:59
  • $\begingroup$ "…only needs to input the vector formed by its derivative into ode45 to solve it" the input of ode45 function is the RHS of the ODE system in standard form, read this post for more info: mathematica.stackexchange.com/a/158519/1871 $\endgroup$
    – xzczd
    Dec 4, 2023 at 2:10
  • 1
    $\begingroup$ @xzczd it looks like OP is using old undocumented interface to ode45 based on this answer at Matlab web site mathworks.com/matlabcentral/answers/…. I see now the wave1D is function name in the code given but it should be defined before using it. When I run their code, I run it from top to bottom one line at time, that is why I got the error from Matlab about udefined wave1D $\endgroup$
    – Nasser
    Dec 4, 2023 at 2:12
  • $\begingroup$ @Nasser If I read it right, OP has contained the definition for wave1D in his post but doesn't format it properly. Just tried editing it a bit. Hopefully I haven't made it worse. $\endgroup$
    – xzczd
    Dec 4, 2023 at 2:18
  • 1
    $\begingroup$ Maybe you mean uflt = Table[uf[i][t], {i, 1, N0}]; vflt = Table[vf[i][t], {i, 1, N0}]; dufl = D[uflt, t]; dvfl = D[vflt, t]; eqs1 = Thread[Equal[dufl, vflt]]; eqs2 = Thread[Equal[dvfl, -a^2 k^2 uflt]];? $\endgroup$
    – Goofy
    Dec 4, 2023 at 3:09

2 Answers 2

10
$\begingroup$

but I can't get the correct result…

I'd argue it's mainly caused by careless mistake… Anyway, we don't have much post solving PDE via FFT (Is this the only one?), so let me do a favor.

2 issues here:

  1. You haven't correctly translated those uts and vts in original MATLAB code. There actually exist 2 different types of ut&vt in the code: one is for storing initial conditions, one is the independent variable of function wave1D. Since you've coded the system as an explicit list of symbolic equation in Mathematica, they need to be identified.

  2. After the FFT via Fourier, we need to inverse it. OP simply doesn't code this part.

The following is a fix with minimal modification for the original code:

L = 80;
N0 = 256;
x = L/N0*Range[-N0/2, N0/2 - 1];
k = (2*Pi/L)*Join[Range[0, N0/2 - 1], Range[-N0/2, -1]];

u = 2*Sech[x];
ut = Fourier[u, FourierParameters -> {-1, 1}];
vt = ConstantArray[0, N0];
uvt = Join[ut, vt];

a = 1;

ufl0 = Table[uf[i][0], {i, 1, N0}];
vfl0 = Table[vf[i][0], {i, 1, N0}];

ufl = Table[uf[i][t], {i, 1, N0}];
vfl = Table[vf[i][t], {i, 1, N0}];

dufl = Table[D[uf[i][t], t], {i, 1, N0}];
dvfl = Table[D[vf[i][t], t], {i, 1, N0}];
eqs1 = Thread[Equal[dufl, vfl]];
eqs2 = Thread[Equal[dvfl, -a^2  k^2  ufl]];
int1 = Thread[Equal[ufl0, ut]];
int2 = Thread[Equal[vfl0, vt]];
rs = NDSolveValue[{eqs1, eqs2, int1, int2}, {ufl, vfl}, {t, 0, 20}];

mat = Table[rs[[1]], {t, 0, 20, 0.5}];

ListLinePlot3D[Chop@InverseFourier[#, FourierParameters -> {-1, 1}] & /@ mat, 
 PlotRange -> All, DataRange -> {{-L/2, L/2}, {0, 20}}]

enter image description here

The code above can of course be simplified, but I'm too lazy to continue. Let me end this answer with a fully NDSolve-based solution for the problem as comparison:

sys = With[{u = u[t, x]}, 
           {D[u, t, t] == a^2  D[u, x, x], 
           {u == 2 Sech[x], D[u, t] == 0} /. t -> 0, 
            (u /. x -> -L/2) == (u /. x -> L/2)}]

sol = NDSolveValue[sys, u, {t, 0, 20}, {x, -L/2, L/2}]

Plot3D[sol[t, x], {t, 0, 20}, {x, -L/2, L/2}, PlotRange -> All, PlotPoints -> 100]

enter image description here


Update

…If you insist on solving the problem based on the same logic as that of MATLAB, then we need a solver similar to ode45 i.e. an ODE solver taking the right hand side of ODE in standard form as the input. Luckily I've already written one here. Let's use it:

L = 80;
N0 = 256;
x = (L Range[-(N0/2), N0/2 - 1])/N0;
k = (2  π )/L Join[Range[0, N0/2 - 1], Range[-(N0/2), -1]];
u = 2 Sech[x];
ut = Fourier[u];
vt = ConstantArray[0, N0];
uvt = Join[ut, vt];
a = 1;
t0 = 0;
tend = 20;    
dt = 1/20;
t = Range[t0, tend, dt];

wave1D = With[{n = N0, a = a, k = k}, 
   Function[{t, uvt}, 
    Module[{ut, vt}, {ut, vt} = Partition[uvt, n]; Join[vt, -a^2  k^2  ut]]]];

(* Definition of rk4stepfunc is not included in this post,
   please find it in the link above. *)
step = (rk4stepfunc /. (Verbatim[_Real]) :> (_Complex))@wave1D;

solmat = FoldList[step[#2[[1]], #, #2[[2]]] &, uvt, 
    Transpose@{Most@t, Differences@t}]; // AbsoluteTiming

ListDensityPlot[Chop[InverseFourier /@ solmat[[;; ;; 10, 1 ;; N0]]], 
 DataRange -> {{-L/2, L/2}, {t0, tend}}, PlotRange -> All]

enter image description here


Update 2

If what you want is to convert the wave1D to Mathematica code and use it in NDSolve, define a black-box function. I'll show a somewhat advanced way for defining it just for fun:

L = 80;
N0 = 256;
x = (L Range[-N0/2, N0/2 - 1])/N0;
k = 2 π/L Join[Range[0, N0/2 - 1], Range[-N0/2, -1]];
u = 2 Sech[x];
ut = Fourier[u];
vt = ConstantArray[0, N0];
uvt = Join[ut, vt];
a = 1;
t0 = 0;
tend = 20;

wave1Dcompiled = 
 With[{n = N0, a = a, k = k}, 
  Compile[{t, {uvt, _Complex, 1}}, 
   Module[{ut, vt}, {ut, vt} = Partition[uvt, n]; Join[vt, -a^2  k^2  ut]], 
   RuntimeOptions -> EvaluateSymbolically -> False]]

solmat = NDSolveValue[{vec'[t] == wave1Dcompiled[t, vec[t]], vec[0] == uvt}, 
    vec /@ Range[0, tend, 1/2], {t, 0, tend}]; // AbsoluteTiming

ListPointPlot3D[Chop[InverseFourier /@ solmat[[All, ;; N0]]], PlotRange -> All, 
 DataRange -> {{-L/2, L/2}, {t0, tend}}]

enter image description here

$\endgroup$
2
  • $\begingroup$ Thank you very much for your answer. Because I tried to solve PDE by FFT for the first time, I didn't know whether the constructed symbolic function was correct or not, so I lost confidence. This question has been bothering me for several days, and I can finally sleep well tonight. Thanks again. $\endgroup$
    – mozeq
    Dec 4, 2023 at 3:50
  • $\begingroup$ Wow, the update is very useful to understand the matlab solve progress. $\endgroup$
    – mozeq
    Dec 4, 2023 at 7:21
2
$\begingroup$

I want to give a more concise code form by learning other people's answers. The original differential equation is not reduced in order.

Definition of fftshift is from the following post:

What's the correct way to shift zero frequency to the center of a Fourier Transform?

L = 80;
N0 = 256;
x = L/N0*Range[-N0/2, N0/2 - 1];
k = (2*Pi/L)*fftshift@Range[-N0/2, N0/2 - 1];
u = 2*Sech[x];
ut = Fourier[u, FourierParameters -> {-1, 1}];
a = 1;
U[t_] := Table[uf[i][t], {i, 1, N0}];
int1 = U[0] == ut;
int2 = U'[0] == ConstantArray[0, N0];
eqs = U''[t] == -a^2 k^2 U[t];
 
rs = NDSolveValue[{eqs, int1, int2}, U[t], {t, 0, 20}];
mat = Table[rs, {t, 0, 20, 0.5}]
ListLinePlot3D[
 Chop@InverseFourier[#, FourierParameters -> {-1, 1}] & /@ mat, 
 PlotRange -> All, DataRange -> {{-L/2, L/2}, {0, 20}}]
$\endgroup$
2
  • $\begingroup$ Adjustions for FourierParameters is actually redundant. $\endgroup$
    – xzczd
    Dec 4, 2023 at 5:55
  • $\begingroup$ Yes, the different FourierParamaters have the same result after InverFourier.But if you want the same result as the fft function in Matlab, you may need this parameter. $\endgroup$
    – mozeq
    Dec 4, 2023 at 7:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.