11
I'm not sure how to explain this behavior, but I've found a solution. Just move the condition to 2nd argument of WhenEvent:
WhenEvent[Mod[t, 2 π] == 0, If[t > 50, Sow@{t, x[t], x'[t]}]]
9
args = Array[x, 20];
D[f @@ args, {args, 1}]
D[f @@ args, {args, 2}]
7
If you write
s = 1; beta = 0.1; M = 0.2; Pr = 1; xmax = 7; (* note xmax, not max *)
then
pfun =
ParametricNDSolveValue[
{f'''[x] + f[x]*f''[x] + beta*(1 - f'[x]^2) + M^2*(1 - f'[x]) == 0,
1/Pr*T''[x] + f[x]*T'[x] == 0,
f[0] == s, f'[0] == L1, T[0] == 1, f'[xmax] == 1, T[xmax] == 0},
{f, T}, {x, 0, xmax}, {L1}];
will return ...
5
This ODE system can be solved using the option Method -> "StiffnessSwitching", although the computation proceeds slowly.
sol = NDSolveValue[{f'''[y] + f[y] f''[y] + 4 - (f'[y])^2 == 0,
g''[y] + 0.01*f[y] g'[y] == 0, f[0] == 0, f'[0] == 0, g'[0] == -1,
f'[10] == 2, g[10] == 0}, {f[y], g[y]}, {y, 0, 10},
Method -> "StiffnessSwitching", ...
4
You could refine the mesh:
test = ParametricNDSolve[{u''[x] - (1/nu)*u[x]*u'[x] == 0,
u[-1] == 1 + delta, u[1] == -1}, u, {x, -1, 1}, {nu, delta},
Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}];
u[0.05, 0.1][0.5] /. test
1.09923
4
As pointed out in the comment, OP's MATLAB implementation doesn't seem to be correct, so I'll just give a solution to OP's problem.
Since the heat flux continuity involves in, there's indeed something new in this problem compared to previous ones. We first handle the continuity condition. Notice that, strictly speaking, 1D heat conduction equation is not
$...
4
OK, write my comment as answer. If I understood OP's intention correctly, he wants to solve for y[x] and plot y'[x]. The equation contains a parameter L, so we need to use ParametricNDSolve:
a = 2.7*10^-3; θ = Pi/6; g = 9.8; ρ = 1000; σ = 70*10^-3;
δp[L_] := 2*σ/L - g*(ρ*L)/2;
sol = ParametricNDSolve[{1 + y'[x]/Sqrt[(1 + y'[x]^2)] ==
x^2/(2*a^2) + (δp[...
3
Quite a few problems here.
your syntax for NDSolve appears incorrect in more ways than one: you will want to look it up in the documentation.
you should use == (double equal) for equations, not = (which is used to set values instead); you will also need to ClearAll[x,y] before re-running your code;
your boundary conditions do not make sense: is e.g. x[t] ...
3
VectorPlot shows you the phase curves of your ode:
vplot = VectorPlot[{1, Sqrt[RealAbs[γ]]}, {t,0, 2}, {γ, -1, 1} ]
Alternativly you might solve the ode numerically
gamma =
ParametricNDSolveValue[
{γ'[t] == Sqrt[RealAbs[γ[t]]], γ[0] == - c^2}, γ, {t,0, 3}, c]
Show[{Plot[Table[gamma[c][t], {c, 0, 1, .1}],{t, 0, 3}],
Plot[ gamma[1][t] , {t, 0, 3},...
2
Try
M = Array[Subscript[y, #1, #2][t] &, {2, 2}];
M0 = {{1, 2}, {-6, 4}};
ci = Thread[Flatten[M] == Flatten[M0]] /. {t -> 0};
DSolve[{i D[M, t] == M/20, ci}, Variables[M], t]
2
We can simplify the procedure since the f-equation can be solved totaly independent of the g-function. Thanks to the inspiration by @bbgodfrey.
First solve the f-equation with initial conditions. (I found parameters i1 and i2 with FindRoot since shooting method did not work with my version 8.0 in that case.)
fsol[i_?NumericQ] :=
f /. First@
NDSolve[{...
2
I cheated a little here.
Btw, I do not think your formula is correct. It should be
(page 167, Differential equations and their applications, 4th edition. By Braun).
Hence, you can do
ClearAll[x,y,k];
sol=y[x]/.First@DSolve[y''[x]+k*y[x]==0,y[x],x];
convert[sol,k,x]
Where
convert[sol_, k_, x_] := Module[{a, b, expr},
expr = TrigExpand[sol];
b = ...
2
You can regard the trz separetly.
First I rationalize all parameters and rewrite equation for simple testing later.
ode = {0 == -τrz[
r] + η/(1 + (ϵ λ/η) τzz) u'[r],
1/r D[r τrz[r], r] == dpdz, u[Ri] == 0, u[Ro] == 0}
Since trz does not depend on u, treat it separatly.
dsol = DSolve[τrz[r]/r + Derivative[1][τrz][r] ==
1000, τrz, r]
{* {{...
1
Alternative numerical formulations as a matrix DE:
ic = RandomReal[1, {2, 2}]; (* set initial condition to whatever is desired *)
Block[{i = 1}, (* set the value of i *)
sol = NDSolve[{i y'[t] == y[t]/20, y[0] == ic}, y, {t, 0, 1}]
]
Block[{i = I}, (* in case i = Sqrt[-1] *)
sol = NDSolve[{i y'[t] == y[t]/20, y[0] == ic}, y, ...
1
Hint.
Solve it as an initial value problem, and then use a shooting method. The following formulation
g = 15/10000;
Ri = 38/1000;
Ro = Ri + g;
m = 1/1000;
ηpos = 1/10;
τy = 2;
ϵ = 1/10;
λ = 1;
dpdz = 10^3;
η = τy/(u'[r] + 10^-10) (1 - Exp[-m u'[r]]) +ηpos;
τzz = τrz[r]^2 (2 λ)/η;
ode = {τrz[r] == η/(1 + (ϵ λ/η) τzz) u'[r], 1/r D[r τrz[r], r] == dpdz,
...
1
There are typos that occurred during the non-dimensioning. If the length of the river is L, then the dimensionless variable x varies from 0 to 1. With this in mind, the code works, but you need to know the boundary conditions.
A = 2100;
v = 4.32*10^4;
k1 = 8.27;
k2 = 44.1;
k = 0.007;
alpha = 16.5;
S = 0.01;
dp = 3.456*10^6;
dx = 3.456*10^6;
L = 325000;
q = ...
1
Here's a way to get the desired plot from NDSolve. An event when g[t] == 0 is reached to make sure the zero solution is followed, followed by an event when the time delay $c^2$ has elapsed:
ClearAll[bvpValue];
bvpValue[g0_?NumericQ, tt_?(VectorQ[#, NumericQ] &), delay_] :=
Module[{t0, t1, t2, ics, ode1, deqns, backsub, start = Infinity},
t0 = ...
Only top voted, non community-wiki answers of a minimum length are eligible
