Does anybody know how to calculate self-calling differential operators fast?

1 view (last 30 days)
Mohammad Shojaei Arani
Mohammad Shojaei Arani on 10 Nov 2021
Commented: Walter Roberson on 10 Nov 2021
Hello friends,
I have a code with a for-loop where inside the for-loop there is a differential operator which uses a function, say A0, as input and performs
some operations and fthe results of the calculations are asigned to A0 again to be used for the next iteration. Intuitively, we understand that
the expressions in the for-loop get more and more complex at each iteration so that the computation time increases.
Long story short: In bellow I attach my code. The double for-loop (with tic and toc) becoms computationlly expensive
very when I increase J and K. I hope you have an idea to imrpve the speed of this code!
Thanks lot in advance!
Babak
clc;
J=4;K=4;
H=sym('H',[J 1]);dt=0.1;
syms x par1 par2 par3 par4 par5
for j=1:J
H(j)=horner(simplify(exp(x^2/2)*diff(exp(-x^2/2),j)),x);
end
syms y y0 positive
muY(y)=-(y*(y*(y*(par1*y*par5^2 - par1*par2*par5) + par2*par3 + par1*par5^2) - par1*par2*par5))/(par2*par5*y^2 + par2*par5);
Hz=subs(H,x,(y-y0)/dt^(1/2));
eta=sym('eta',[J 1]);
tic;
for j=1:J
A0=Hz(j);S=A0;
for k=1:K
A0=muY(y)*diff(A0,y)+1/2*diff(A0,y,2);
S=S+A0*dt^k/factorial(k);
end
eta(j)=1/factorial(j)*limit(S,y,y0);
end
toc;
  1 Comment
Walter Roberson
Walter Roberson on 10 Nov 2021
I wonder...
If you were to use
syms MUY(y)
and
A0=MUY(y)*diff(A0,y)+1/2*diff(A0,y,2);
then you should be able to get out a response in terms of powers of MUY(y) and derivatives of MUY(y) . At that point, you can use findSymType() to locate derivatives of MUY, and find the maximum derivative, as well as building up a list of which derivatives are used. Loop through iteratively building up the derivatives to the maximum used. Then subs() for the derivative actually used in S.
I think doing this should lower the cost of doing the derivatives of muY. The subs() at the end might be expensive, but I suspect you will still end up saving time.

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by