non-uniform grid used in a for loop

Hi,

I am trying to integrate an ODE (an IVP, dR/dS = ..., R(0)=0) via trapezoidal rule, and the way I implemented it is via a for loop j = 1:1:J, and J is the number of points I use to discretize the range of S (0 to 200) over which the integration takes place. I first used a uniform discretization for S, and the code take the rough form of below:

dS = (S_max-S_min)/(J-1) % J points hence (J-1) intervals

R(1) = 0; A(1) = 0; B(1) = 0; % initial value of R, A and B

for j = 2:1:J

A(j) = (2/(vm*(j-1)*dS));

B(j) = (2*(j-1)*dS)/(vm*((j-1)*dS)^2);

a = A(j)/2;

b = 1/dS + B(j)/2;

c = A(j-1)*R(j-1)/2 - (1/dS - B(j-1)/2)*R(j-1) - 1;

R(j) = -c/(b + sqrt(b^2-4*a*c));

end

I believe that applying a non-uniform discretization for the range of S would render better results, more specifically to have more points around S = 100. A very manual and probably dumb way that I thought of was to predefine the number of points I want when its around S = 100 and when it's not. for example:

dS_1 = (90-S_min)/100; % 100 points from S_min to 90

dS_2 = (110-90)/200; % 200 points from 90 to 110

dS_3 = (S_max - 110)/100; % 100 points from 110 to S_max

Once the new dS values are defined, I used them in 3 different for loops (Smin to 90, 90 to 110, 110 to Smax) like this

R(1) = 0; A(1) = 0; B(1) = 0; % initial value of R, A and B

for j = 2:1:100 % the first 100 points

A(j) = (2/(vm*(j-1)*dS_1));

B(j) = (2*(j-1)*dS_1)/(vm*((j-1)*dS_1)^2);

a = A(j)/2;

b = 1/dS_1 + B(j)/2;

c = A(j-1)*R(j-1)/2 - (1/dS_1 - B(j-1)/2)*R(j-1) - 1;

R(j) = -c/(b + sqrt(b^2-4*a*c));

end

I then repeat this kind of method for j = 101:1:300 and j = 301:1:400, using dS_2 and dS_3 respectively.

My question is that, when I ran the code, I found that it is very inefficient and slow (the codes above are simplified version, the actual one is much longer and complex involving more parameters), so could someone maybe share with me some other possible ways I may do some research on to implement this non-uniform grid for this integration process?

Thank you for your time!

[Torsten]

Why don't you use one of the integrators from the Matlab solver suite ? It adapts the step size automatically in order to ensure a prescribed error tolerance.

Look up ode45, ode15s, e.g.

[Mingze Yin]

Hi, yes that's what I used at first, but in the context of my research topic, there are some other issues in the larger setting if I directly use ode solvers. The trapezoidal rule is the more sensible method in my case.

[Torsten]

In my opinion, it does not make much sense to predefine a non-uniform grid. If model parameters change, it might be necessary to alter the values 90 and 110. The usual procedure is to adapt step size during the integration process. This is more efficient than working with a uniform stepsize over the complete integration interval and lowers the error in solving the equations. There should be adaptive integrators for the trapezoidal rule in the internet. You should give it a try.

[Mingze Yin]

oh ok, thanks for the info, I'll look into it!

[Mingze Yin]

Hi, I did some research on the adaptive trapezoidal rule and what I understand is that there is always the comparison between the estimated value and the "real" value?(correct me if I'm wrong) However, what if I do not have a reference value to compare to in the first place? Is it still possible to create the adaptive grid?

[Torsten]

You don't compare with the real value, but you compare the results obtained if you compute u(t+deltat) once with stepsize deltat and once with two times step size deltat/2. If the two values you receive are close, you can choose a larger stepsize for the next step. If they differ much, you have to reduce the step size.

But you shouldn't try to implement adaptive stepsize control on your own. You should use a reliable solver instead.

This might interest you:

sciencedirect.com/science/article/pii/037704279400007N