How can the pdepe solver solution be used as a new initial condition for the next iteration? i.e using pdepe within a loop

I solved the convection diffusion equation using pdepe solver.

I would like to solve this equation within a loop. i.e

the output solution at the tout is a vector, which should be used as new initial condition.

It would be appreciated to help in this issue.

[maria atiq]

Hey, were you able to solve this issue? I want to use my pdepe solver within a loop where i have designed a PID to control the output of my pde. Kindly share if you were able to run the pdepe as a loop.

Answer 1

[Venkatachala Sarma]

Hi Younus,

This could be done by using a variable shared between the function calling 'pdepe' and the function pdex1ic(x) which gets called to manage the initial conditions. Let us take a global variable as an example.

Consider the first example in the following documentation page.

http://in.mathworks.com/help/matlab/ref/pdepe.html

In this example, you could declare a global variable which will be updated with the solution of 'pdepe' every time it runs. Just equate this global variable to the output of the function pdex1ic which manages the initial condition.

I am pasting the modified version of pdex1 and pdex1ic functions alone. Remaining of the example code is retained.

if true
  function pdex1
    global takeOut;
    m = 0;
    x = linspace(0,1,20);
    t = linspace(0,2,5);
    takeOut = 0;
    for i = 1:5
       sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
       u = sol(:,:,1);
       takeOut = u(end,end);   
       surf(x,t,u) 
       title('Numerical solution computed with 20 mesh points.')
       xlabel('Distance x')
       ylabel('Time t')
       figure
       plot(x,u(end,:))
       title('Solution at t = 2')
       xlabel('Distance x')
       ylabel('u(x,2)')
    end
  end
end

The above function puts the pdepe function in a loop and uses the global variable named takeOut to carry the current solution. This is just to show an example and you can explore the other possibilities. The pdex1ic function shown below just copies the value in this variable.

if true
  function u0 = pdex1ic(x)
    global takeOut;
    u0 = takeOut;
  end
end

Update these functions in the example and run it. You could see that the outputs vary depending upon the initial conditions.

Hope this helps.

Regards Venkatachala Sarma

[ANDREAS CHARALAMBOUS]

Dear Sarma,i have a very similar problem Younus and your help will be much aprreciated.I am trying to simulate a process with five steps and i have 7 variables which calculated with pdepe(u1,u2,u3,u4,u5,u6,u7) and i want the solution profiles u1(end,:),u2(end,:),u3(end,:),u4(end,:),u5(end,:),u6(end,:),u7(end,:) to be my initial conditions for the next step.So for example the solution profiles from the first step will be the initial conditions for the second step.The solution profiles for the second step will be the initial conditions for the third step....and the solution profiles from the fifth step will be the initial conditions for the first step.So i want to create a cycle like this.How can i implementing this in pdepe?And the second problem is that for each step i have different boundary conditions.

[Antonio López]

I have exactly the same problem. I need u1(end,:) to be the initial condition of next PDE. Obviously, u1(end,:) is a vector of same length as x (spatial coordinate).

In the solution that Sarma gave to Younus, takeOut is only a number, that's why it works without any problem.

But Andreas and I need the whole profile u1(end,:), not just a number. In other words, we need to asign a different value to every x position (discrete) in our initial condition.

[Antonio López]

Dear Andreas,

I've found a solution.

Let me denote x as the spatial coordinate. For example, imagine that x is disretised between 0 and 0.05 with intervals of 0.01, so x=0:0.01:0.05. The lenght of this vector is 6.

From the solution of our first PDE, we obtain u1. For example, it could be that u1(end,:)=[1,2,3,4,3,1] (I'm inventing the numbers). Of course, the length of u1(end,:) is 6.

First of all, you have to copy your spatial coordinate discretization to another variable, for example, xcopy. This variable must be declared as global so that it could be accesible from the initial condition function.

In the initial condition function of the next PDE, you can do this:

function u0 = pdexAAAAA_ic(x)

       global xcopy;
       u0 = u1(end, find(xcopy==x,1,'first'));

end

The goal of this initial condition function is just to fill the first row of the next PDE solution, so internally it's like a for loop in which x is the index of the loop. Every time u0 is called, you are finding the position of the current value of x in the vector xcopy. Doing this, you can access the previous solution u1 in the time end and in the correct position.

I hope it would be helpful. Inserting a breakpoint in the initial condition function would help to understand the process during the execution.

[Haonian Shu]

Dear Antonio,

I am have the same issue and have been trying to follow your solution. Yet I cannot figure out where should I declare the solution u1, since u1 is not global, how can it be extracted by the ic function. also, I am struggling to understand the goal of this initail condition function, the goal of xcopy. actually, all i want to find is to give pdepe a custumized, discrete (arbitrary) initial condition rather than a number or a simple analytical function.

Best