## 195 Reputation

10 years, 25 days

Amir

## reverse solution of diffusion equation...

Maple 2019

Hello all

I am going to find the diffusion constant in the following equation :

PDE:=diff(C(x,t),t)=d*diff(C(x,t),x,x);

In a way that ,

evalf(Int(C(x,t=specific),x=0..L))/L - m_number(t=specific)=0

In other words, It would be an iterative procedure to guess "d" and then find the solution C (something like shooting method). The objective is to find "d" in a way that the average of C (solution of PDE) would be equal to m_number at that specific time.

Also, Since I have 5 specific time and m_number, I will have 5 different d. So, I need to use least square method or other optimization technique find one finalized d value.

inverse_numerical_diffusion.mw

I am not good at using Proc command and  I think the error of my code is because of that.

Would you help me to find out my problem in the code and any new idea for solviong this equaltion would be greatly appreciated.

Amir

## Diffusion Equation - inverse solution fo...

Maple 2019

Hello everyone

I have the solution of diffusion equation from Help of maple website. I put the code here

*****************************

restart: with(plots):

unprotect(D);

alias(c[0]=c0, c[1]=c1, c[2]=c2);
PDE:=diff(C(x,t),t)=D*diff(C(x,t),x,x);
IBC:={C(x,0)=cx0, C(0,t)=ct0, D[1](C)(10,t)=0};
ct0:=1;
cx0:=0;
D:=1;
pds:=pdsolve(PDE,IBC,numeric);
L1:=[0.01, 0.1, 1, 5, 10];
L2:=[red, green, yellow, blue, magenta, black];
for i from 1 to 5 do
pn[i] := pds:-plot(t=L1[i], color=L2[i]):
end do:
display({seq(pn[i], i=1..5)}, title=`Numerical solution at t=0.01, 0.1, 1, 5, 10`);

****************************

the code is working perfectly. But, My question is how can I found the diffusion constant (D) if I have the solution ( C(x,t) ).  Probably it should be an algorithm which use least square method to find (D) based on the data C(x,t).

I am looking for a fast and efficient algorithm if there is any.

Sincerely yours,

Amir

## differentiate in Series...

Maple 18

Hi everyone

I have a problem with differentiating an equation. I have the following code:

(sum((x(i)-M)^3, i = 1 .. n))/(n*((sum((x(i)-M)^2, i = 1 .. n))/(n-1))^(3/2));

The values of x(i) is available in Excel and I want to import that. For this case, I copy the column and paste in maple and it seems work fine. But, I want to calculate the differentiate of the summation with respect to x(i) and M and find the answer. I don't know how to write the code appropriately.

In other words, how could I differetiate the following function with respect to x(i) and put the values from excel into it.

(sum((x(i)-M)^3, i = 1 .. n))/(n*((sum((x(i)-M)^2, i = 1 .. n))/(n-1))^(3/2));

## Convergence problem of ODE system...

Maple

Dear collegues

Hope you are fine

I wrote a code to solve a system of ODEs.

The code solve the problem for higher values of parameter NBT>=5. When I decrease it to NBT=0.2, the code didnt converge. I did my best but I couldnt get the results.

I would be most grateful if you help me at this problem

The code is attached

Thank you

Final_code.mw

Amir

## Convergence problem for some values of p...

Maple

Dear collagues

Hi,

I write a code to solve a system of ODE. It solve the ODES in a wide range of parameters but as I decrease NBT below 0.5, it doesnt converge. I do my best but I couldn't find the answer. Would you please help me? Thank you

Here is my code and it should be run for 0.1<NBT<10. the value of NBT is input directly in res1.

restart:
EPSILONE:=1000:
Digits:=15:

a[mu1]:=5.45:
b[mu1]:=108.2:
a[k1]:=1.292:
b[k1]:=-11.99:

rhop:=4175:
rhobf:=998.2:
mu1[bf]:=9.93/10000:
k1[bf]:=0.597:

rhost(eta):=1-phi(eta)+phi(eta)*rhop/rhobf;
k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta);

eq1:=(diff(u(eta), eta))*a[mu1]*(diff(phi(eta), eta))+2*(diff(u(eta), eta))*b[mu1]*phi(eta)*(diff(phi(eta), eta))+((diff(u(eta), eta))+(diff(u(eta), eta))*a[mu1]*phi(eta)+(diff(u(eta), eta))*b[mu1]*phi(eta)^2)/(eta+EPSILONE)+diff(u(eta), eta, eta)+(diff(u(eta), eta, eta))*a[mu1]*phi(eta)+(diff(u(eta), eta, eta))*b[mu1]*phi(eta)^2+1-phi(eta)+phi(eta)*rhop/rhobf:
eq2:=(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2)*(diff(T(eta), eta))/(eta+EPSILONE) + (a[k1]*(diff(phi(eta), eta))+2*b[k1]*phi(eta)*(diff(phi(eta), eta)))*(diff(T(eta), eta))+(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2)*(diff(T(eta), eta, eta));
eq3:=diff(phi(eta),eta)+phi(eta)/(N_bt*(1+gama1*T(eta))^2)*diff(T(eta),eta):

eq2:=subs(phi(0)=phi0,eq2):
eq3:=subs(phi(0)=phi0,eq3):

eval(dsolve({eq3,phi(1)=phiv},phi(eta)),(T)(1)=1):
Phi:=normal(combine(%)):
Teq:=isolate(eval(eq2,Phi),diff(T(eta),eta)):
ueq1:=eval(eq1,Phi)=0:
ueq2:=subs(Teq,ueq1):

lambda:=0;
Ha:=0;
N_bt:=cc*NBT+(1-cc)*0.8;
kratio:=k1[p]/k1[bf]:

GUESS:=[T(eta) =0.0001*eta, u(eta) =0.1*eta, phi(eta) = 0.3*(eta-1)^4];
res1 := dsolve(subs(NBT=0.48,gama1=0.2,phiv=0.06,{eq1,eq2,eq3,u(0)=lambda*D(u)(0),D(u)(1)=0,T(0)=0,phi(1)=phiv,T(1)=1}), numeric,method=bvp[midrich],maxmesh=4000,approxsoln=GUESS, output=listprocedure,continuation=cc):
G0,G1,G2:=op(subs(subs(res1),[phi(eta),u(eta),diff(T(eta),eta)])):

masst:=evalf(int((1-G0(eta)+G0(eta)*rhop/rhobf)*G1(eta), eta = 0..1));
heatt:=(1+a[k1]*G0(0)+b[k1]*G0(0)^2)*G2(0);

plots:-odeplot(res1,[[eta,T(eta)]],0..1,legend=[T],color=["Black"],linestyle=Solid,axes=boxed,thickness=3);
plots:-odeplot(res1,[[eta,u(eta)]],0..1,legend=[u],color=["Black"],linestyle=Solid,axes=boxed,thickness=3);
plots:-odeplot(res1,[[eta,phi(eta)]],0..1,legend=[phi],color=["Black"],linestyle=Solid,axes=boxed,thickness=3);

>
>

Thank you

Amir

 1 2 3 4 5 6 7 Last Page 1 of 11
﻿