# How to solve numerically system of linear coupled PDE using Finite Difference? [closed]

How to solve this coupled linear equation with the 2D plot of the variation of parameter with constant with I.C: u(y,0)=c(y,0)=theta(y,0)=0, B.C: u(0,t)=1,c(0,t)=theta(0,t)=t.

• what values do the other parameters Gr,Gm,R,Pr,Sc have? Oct 26, 2022 at 20:36
• And, is this a question about the software Mathematica? Oct 27, 2022 at 2:33
• These equations may be solvable symbolically, because they are linear and autonomous. Oct 27, 2022 at 4:56
• How ? can you explain the code, please? Oct 27, 2022 at 5:03
• @AshirbadRath please stop deleting content of this question, it makes the topic useless.
– Kuba
Oct 27, 2022 at 13:19

Since this forum is for using Mathematica to solve problems, here is how you would solve these using NDSolve.

If you want code using finite difference, I suggest the forum https://scicomp.stackexchange.com/ instead.

I fixed few things. First your BC and IC is not consistent, You say u(y,0)=0 but then say u(0,t)=1. So I changed u(0,t)=0.

You also do not have enough spatial B.C. You have derivatives w.r.t. y twice, but only provided one BC. So I added more BC's to make NDSolve happy. You can change these as needed.

You also did not provide numerical values for Gr,Gm,RmPr,Sc do I added dummy values for these. You can change them to the actual values.

ClearAll[u, y, t, c, θ, Gr, Gm, R, Pr, Sc]
pde1 = D[u[y, t], t] == Gr* θ[y, t] + Gm*c[y, t] + D[u[y, t], {y, 2}]
pde2 = D[θ[y, t], t] == 1/Pr*D[ θ[y, t], {y, 2}] - R/Pr*θ[y, t];
pde3 = D[c[y, t], t] == 1/Sc*D[c[y, t], {y, 2}]

ic = {u[y, 0] == 0, θ[y, 0] == 0, c[y, 0] == 0}

bc = {u[0, t] == 0, θ[0, t] == t, c[0, t] == t,
(D[u[y, t], y] /. y -> 0) == 0, (D[θ[y, t], y] /. y -> 0) == 0,
(D[c[y, t], y] /. y -> 0) == 0};

parm = {Gr -> 1, Gm -> 2, R -> 3, Pr -> 4, Sc -> 5};

NDSolve[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ, c}, {t, 0, 2}, {y, 0, 1}];


Update to make plot

You can make plots as follows

{usol, θsol, Csol} =
NDSolveValue[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ,
c}, {t, 0, 2}, {y, 0, 1}]
Plot3D[usol[t, y], {t, 0, 2}, {y, 0, 1}, PlotLabel -> "u solution",
AxesLabel -> {"time", "y", "u(t,y)"}, BaseStyle -> 14]


Plot3D[θsol[t, y], {t, 0, 2}, {y, 0, 1},
PlotLabel -> "θ solution",
AxesLabel -> {"time", "y", "θ(t,y)"}, BaseStyle -> 14]


and so on. To plot one solution against the other, you could look into ParametricPlot and ParametricPlot3D. I am not sure what you want to plot, but there are examples how to use these commands in help.

Update to plot with varying time

like when varying t other parameters are constant.

You could use Manipulate for this. Here is an example for u(y,t) which you can the code for the other solutions. To play with varying the other parameters, you just make a slider for each one and then solve again whenever any parameter changes.

Manipulate[
Grid[{
{Row[{"Time = ", t0, " seconds"}]},
{Plot[usol[y, t0], {y, 0, 1},
AxesLabel -> {"y", "u(t,y)"}, BaseStyle -> 12,
ImageSize -> 300, PlotRange -> {Automatic, {-.1, .3}},
PlotStyle -> Red, GridLines -> Automatic,
GridLinesStyle -> LightGray]
}}],
{{t0, 0, "time"}, 0, 2, .01, Appearance -> "Labeled"},
TrackedSymbols :> {t0}
]


how I get 3 different line accordance to this 3 change value in R

You can make new Manipulate variable for R and solve new each time this is clicked. Like this

Make similar variable for the other variables. Here is the new code

Manipulate[
Module[{u, y, t, c, θ, Gr, Gm, Pr, Sc, pde1, pde2, pde3},
pde1 =
D[u[y, t], t] ==
Gr*θ[y, t] + Gm*c[y, t] + D[u[y, t], {y, 2}];
pde2 =
D[θ[y, t], t] ==
1/Pr*D[θ[y, t], {y, 2}] - R/Pr*θ[y, t];
pde3 = D[c[y, t], t] == 1/Sc*D[c[y, t], {y, 2}];
ic = {u[y, 0] == 0, θ[y, 0] == 0, c[y, 0] == 0};
bc = {u[0, t] == 0, θ[0, t] == t,
c[0, t] == t, (D[u[y, t], y] /. y -> 0) ==
0, (D[θ[y, t], y] /. y -> 0) ==
0, (D[c[y, t], y] /. y -> 0) == 0};
parm = {Gr -> 1, Gm -> 2, Pr -> 4, Sc -> 5};
{usol, θsol, Csol} =
NDSolveValue[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ,
c}, {y, 0, 1}, {t, 0, 2}];
Grid[{
{Row[{"Time = ", t0, " seconds"}]},
{Row[{"R = ", R}]},
{Plot[usol[y, t0], {y, 0, 1}, AxesLabel -> {"y", "u(t,y)"},
BaseStyle -> 12, ImageSize -> 300,
PlotRange -> {Automatic, {-.1, .3}}, PlotStyle -> Red,
GridLines -> Automatic, GridLinesStyle -> LightGray]
}}]
],
{{t0, 0, "time"}, 0, 2, .01, Appearance -> "Labeled"},
{R, {3, 2, 9}},
TrackedSymbols :> {t0, R}
]


How to manipulate simulationsly all three plots for R .

One way is to make 3 plots and then use show. This could be improved but gives you the idea

Updated code

Manipulate[
Module[{u, y, t, c, θ, Gr, Gm, Pr, Sc, pde1, pde2, pde3, R,
parm},
pde1 = D[u[y, t], t] == Gr*θ[y, t] + Gm*c[y, t] + D[u[y, t], {y, 2}];
pde2 = D[θ[y, t], t] == 1/Pr*D[θ[y, t], {y, 2}] - R/Pr*θ[y, t];
pde3 = D[c[y, t], t] == 1/Sc*D[c[y, t], {y, 2}];
ic = {u[y, 0] == 0, θ[y, 0] == 0, c[y, 0] == 0};
bc = {u[0, t] == 0, θ[0, t] == t,
c[0, t] == t, (D[u[y, t], y] /. y -> 0) ==
0, (D[θ[y, t], y] /. y -> 0) ==
0, (D[c[y, t], y] /. y -> 0) == 0};
parm = {Gr -> 1, Gm -> 2, R -> 3, Pr -> 4, Sc -> 5};
{usol1, θsol1, Csol1} = NDSolveValue[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ, c}, {y, 0, 1}, {t, 0, 2}];
parm = {Gr -> 1, Gm -> 2, R -> 30, Pr -> 4, Sc -> 5};
{usol2, θsol2, Csol2} = NDSolveValue[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ, c}, {y, 0, 1}, {t, 0, 2}];
parm = {Gr -> 1, Gm -> 2, R -> 200, Pr -> 4, Sc -> 5};
{usol3, θsol3, Csol3} = NDSolveValue[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ, c}, {y, 0, 1}, {t, 0, 2}];
Grid[{
{Row[{"Time = ", t0, " seconds"}]},
{Plot[{usol1[y, t0], usol2[y, t0], usol3[y, t0]}, {y, 0, 1},
AxesLabel -> {"y", "u(t,y)"}, BaseStyle -> 12, ImageSize -> 300,
PlotRange -> {Automatic, {-.1, .3}},
PlotStyle -> {Red, Blue, Green}, GridLines -> Automatic,
GridLinesStyle -> LightGray,
PlotLegends -> {"R=3", "R=30", "R=200"}]
}}]
],
{{t0, 0, "time"}, 0, 2, .01, Appearance -> "Labeled"},
TrackedSymbols :> {t0}
]

• Thank you so much, up to that I'm clear. But one thing I ask is a 2D plot with various parameter values. Gm has multiple values while others are the constant mean same value. Oct 27, 2022 at 4:05
• Plot[Evaluate[{u[y, t]}, {\[Theta][y, t]}, {c[y, t]} /. sol], {t, 0, 1}, {y, 0, 1}, PlotStyle -> Automatic] . In this plot, I get t an error like Plot::nonopt: Options expected (instead of {y,0,1}) beyond position 2 in Plot[{u[y,t]},{\[Theta][y,t]},{{InterpolatingFunction[{{0.,1.},{0.,2.}},{5,5,1,{25,64},{6,4},0,0,0,0,Automatic,<<3>>},{{0.,<<9>>,<<15>>},{0.,<<9>>,<<54>>}},{DeveloperPackedArrayForm,{0,3,6,9,12,15,18,21,24,27,0072*10^15,1.,0.000033917,<<4790>>}},{Automatic,Automatic}][y,t]}},<<1>>,{y,0,1},PlotStyle->Automatic]. An option must be a rule or a list of rules.` Oct 27, 2022 at 4:30