Mathematica For Loop to Solve PDE: StressIC Function and NDSolve Error Fix

[as144]

I am trying to write a for loop in Mathematica to solve a PDE , But I am facing a problem. My For loop is:
StressIC[s11_, s22_] :=
Limit[(1/((2 Pi)*c*c))*
Exp[-(1/2)*(((s11 - 0)/c)^2 + ((s22 - 0)/c)^2)], c -> 0.001];
countmax = 10;
t =.;
count = 1;
For[count = 1, count <= countmax, count = count + 1, Print[count];
t =.;
ElasticPlasticEq = {D[AC11p[s11, s22]*P[s11, s22, t], s11] +
2*D[AC22p[s11, s22]*P[s11, s22, t], s22] + D[P[s11, s22, t], t] -
D[DC1111p[s11, s22]*P[s11, s22, t]*t, {s11, 2}] -
2*D[DC2222p[s11, s22]*P[s11, s22, t]*t, {s22, 2}] == 0,
P[s11, s22, 0] == StressIC[s11, s22],
AC11p[LB, s22]*
P[LB, s22, t] - (D[DC1111p[s11, s22]*P[s11, s22, t]*t, s11] /.
s11 -> LB) == 0,
AC11p[RB, s22]*
P[RB, s22, t] - (D[DC1111p[s11, s22]*P[s11, s22, t]*t, s11] /.
s11 -> RB) == 0,
2*AC22p[s11, LB]*P[s11, LB, t] -
2*(D[DC2222p[s11, s22]*P[s11, s22, t]*t, s22] /. s22 -> LB) == 0,
2*AC22p[s11, RB]*P[s11, RB, t] -
2*(D[DC2222p[s11, s22]*P[s11, s22, t]*t, s22] /. s22 -> RB) == 0};
ElasticPlasticSol =
NDSolve[ElasticPlasticEq,
P[s11, s22, t], {s11, LB, RB}, {s22, LB, RB}, {t, 0, 0.00000002},
MaxSteps -> {150, 150, 1000000}, StepMonitor :> Print[t]];
t = 0.00000002;
StressIC[s11_, s22_] = Evaluate[P[s11, s22, t] /. ElasticPlasticSol];
Print[StressIC[s11, s22]];
Plot3D[StressIC[s11, s22], {s11, -0.005, 0.005}, {s22, -0.005,
0.005}, AxesLabel -> {"sigma11", "sigma22",
"P[s11,s22,0.00000002]"}, PlotRange -> All] // Print
]

I can have the solution for the first step, but not for the other steps.
the error is:

NDSolve`FiniteDifferenceDerivativeFunction::ddim: Data {<<1>>,<<9>>,<<19>>} is not a rectangular tensor with dimensions {29,15}.

please help me in this matter!

 

[Bill Simpson]

Not having any of the data or the values for a variety of constants in your code or some of your functions makes this difficult to diagnose at best.

I recommend you insert some print statements to show the values for

ElasticPlasticEq, P[s11, s22, t], {s11, LB, RB}, {s22, LB, RB},

just before your ElasticPlasticSol = NDSolve[] line and look carefully at the dimensions of those.

I am guessing that your assignment to ElasticPlasticEq is not what either you or NDSolve is expecting for at least one of your iterations.

If you get that error for every iteration that probably points at a more general problem with your ElasticPlasticEq. If you only get that for some iterations then it is more likely that this fails for some subset of parameters.

I'm guessing if you carefully count the number of elements in each row and column you will the mismatch that generates the error.

 

[as144]

Thank you for your reply, but when I print ElasticPlasticEq even for specific s11, s22 and t, I am getting this:

s11 = 0.01;
s22 = 0.01;
t = 0.00000001;
Print[ElasticPlasticEq];

{0. P[0.01,0.01,1.*10^-8]+(P^(0,0,1))[0.01,0.01,1.*10^-8]+2 (0. P[0.01,0.01,1.*10^-8]+4.88591 (P^(0,1,0))[0.01,0.01,1.*10^-8])-1.25*10^-7 (-68467.2 P[0.01,0.01,1.*10^-8]+0. (P^(0,1,0))[0.01,0.01,1.*10^-8]+0.977181 (P^(0,2,0))[0.01,0.01,1.*10^-8])-9.77181 (P^(1,0,0))[0.01,0.01,1.*10^-8]-2.5*10^-7 (-68467.2 P[0.01,0.01,1.*10^-8]+0. (P^(1,0,0))[0.01,0.01,1.*10^-8]+0.977181 (P^(2,0,0))[0.01,0.01,1.*10^-8])==0,P[0.01,0.01,0]==5.92068*10^-39,-1.49171 P[-0.066,0.01,1.*10^-8]-3.72926*10^-8 (P^(1,0,0))[-0.066,0.01,1.*10^-8]==0,-9.96507*10^-7 P[0.19,0.01,1.*10^-8]-2.49128*10^-14 (P^(1,0,0))[0.19,0.01,1.*10^-8]==0,1.4917 P[0.01,-0.066,1.*10^-8]-2 (5.80742*10^-7 P[0.01,-0.066,1.*10^-8]+9.32316*10^-9 (P^(0,1,0))[0.01,-0.066,1.*10^-8])==0,9.96512*10^-7 P[0.01,0.19,1.*10^-8]-2 (-1.34035*10^-12 P[0.01,0.19,1.*10^-8]+6.2282*10^-15 (P^(0,1,0))[0.01,0.19,1.*10^-8])==0}

how can I follow this to see the dimension?

 

[Dale]

Just count. It is a 6-element array. If the function is expecting a 29 by 15 element matrix then that would be the source of your trouble.

 

[Bill Simpson]

If you look closely at that you also have some very suspicious things.

P[0.01,0.01,1.*10^-8] looks like a function that accepts 3 arguments.

(P^(0,0,1))[0.01,0.01,1.*10^-8] but this looks like I don't know what.

It looks like my suggestion of what you should print and inspect has certainly uncovered at least the first layer of things that have to be fixed.